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Pref ace 


Ify  earliest  memory  at  mathematics  recalls  writing 
numbers  in  the  spaces  of  a  grid.  Little  did  the  first 
grader  know  that  by  creating  this  "box  of  numbers",  he  was 
laying  the  groundwork  for  subsequent  study  in  matrices, 
themselves  very  special  "boxes  of  numbers."  After  studying 
matrix  theory  as  a  high  school  freshman,  I  remember  hoping 
for  a  chance  someday  to  pursue  it  in  great  depth.  Though 
the  hope  has  now  been  realized,  it  will  by  no  means  be  put 
to  rest:  for  my  master's  in  mathematics,  I  will  pick  up 
where  this  thesis  left  off. 

Fortune  doesn't  always  favor  us,  and  as  Machiavelli 
says,  will  strike  after  being  wooed  and  cajoled.  However, 
she  must  have  been  in  a  good  mood  in  the  Fall  of  '84  because 
our  class  had  Dr.  John  Jones  Jr  for  Numerical  Analysis.  His 
presentation  of  matrix  algebraic  concepts  delighted  me,  and 
upon  discovering  his  abiding  interest  and  renown  with  matrix 
equations,  decided  to  study  under  him.  Undoubtedly, 

Dr  Jones  was  one  of  two  teachers  having  a  profound  impact  on 
my  mathematical  development.  1  couldn't  have  wished  for  a 
finer  advisor  and  coach  than  him.  His  patience  and  positive 
outlook  sustained  me  during  those  times  when  I  was  setting 
back  mathematical  science  instead  of  advancing  it.  For  him, 
theory  and  applications  are  one,  though  primacy  goes  to 
theory.  Just  as  it  should  be.  Thanks  for  your  interest  in 
me.  Dr  Jones.  I  will  return  your  concern  years  hence  to  one 


of  my  students. 

To  LTC  (and  Dr.)  Charles  Ebeling,  my  reader  in  the 
Operations  Research  Department,  I  want  to  thank  for  being 
(I  believe)  a  mathematician  at  heart,  and  hence  an  ally  in 
the  department.  If  I  couldn't  have  done  a  mathematical 
thesis,  I  would  have  done  a  simulation  one  and  sought  you 
out  as  an  advisor.  Thanks  for  not  discouraging  me  from 
pursuing  an  interest  begun  years  ago. 
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Abstract 

This  paper  seeks  the  solutions  to  a  system  o-f  equations 
(equalities)  in  n  variables  by  expressing  the  system  in 
matrix  algebraic  -form.  Properties  o-f  the  solutions  to  the 
ensuing  matrix  equation  are  investigated  using  similarity 
transformations.  The  three  types  o-f  matrix  equations  to  be 
studied  are: the  linear  equation - 

AX  =  b 

the  Lypanov  equation 

AX  -  XB  =  C 

the  second-order  Riccati  equation  — 

XDX  +  AX  +  XB  +  C  *  O' 
and  the  third-order  Riccati  equation  — 

XAXBX  +  XCX  +  DX  +  XE  +  F  =  0  , 

The  entries  o-f  all  matrices,  including  the  solution  X,  are 
restricted  to  being  polynomials  in  r  having  complex 
coefficients,  where  r  is  the  n-tuple  of  indetermi nates. 

That  is,  all  matrices  are  elements  of  the  ring  <tmX"(r)  -for 
m  and  n  of  appropriate  size. 

Because  adding  and  multiplying  matrices  (having 
multivariate  polynomial  entries)  is  tedious  in  practice,  an 
interactive  BASIC  program  is  presented  in  the  appendix. 

This  program,  which  can  be  used  on  a  personal  computer, 
permits  the  user  to  perform  operations  on  matrices  having 
multivariate  polynomial  entries.  Via  menu  selections,  the 
user  may  perform 

-weighted  addition  between  two  matrices 


I  Overview  and  Literature  Review 

Overview.  This  paper  seeks  the  solutions  of  the 
third-order  Riccati  matrix  equation 

XAXBX  +  XCX  +  DX  +  XE  +  F  =  0  (1.1) 

where  the  entries  o-f  the  matrices  A,B,C,D,E,F,X  are 
multivariate  polynomials.  In  solving  (1.1)  the  linear 
matrix  equation 

AX  =  b 

the  Lypanov  matrix  equation 

AX  -  XB  *  C 

and  the  second-order  Riccati  matrix  equation 

XDX  +  AX  +  XB  +  C  -  0 

will  be  addressed  because  the  -form  o-f  their  solutions  will 

hint  at  the  nature  o-f  those  to  (1.1).  Because  cursory 

knowledge  of  matrix  algebra  is  sufficient  to  motivate  the 

paper's  thesis — the  solution  to  (1.1) — the  rigorous  (and 

lengthy)  definition  of  terms  and  statement  of  objective  will 

be  made  in  the  next  chapter. 

Why  a  matrix  equation  (such  as  those  given  above)  is 

worthy  of  attention,  let  alone  finding  its  solution,  is  a 

legitimate  concern.  Practical  problems  often  arise  which 

require  the  solution  to  a  system  of  m  equations  in  n 

unknowns,  e.g.,  the  system  of  equations 

ax  +  by  =  c  (1.2) 

dx  +  ey  *  f 


where  all  variables  except  x  and  y  are  known.  One  way  to 


-find  the  solutions  x  and  y  is  to  express  the  system  in  its 
equivalent  matrix  -form 


i:  :i  [3 


(which  is  a  linear  matrix  equation)  and  solve  -for  x  and  y, 


13  -  c  :i  □ 


Thus,  a  system  of  equations  may  give  rise  to  a  matrix 
equation  whose  solution  in  turn  gives  the  answer  to  the 
original  system. 

Other  more  complex  systems  may  have  a  matrix 
representation.  For  instance 


[x  y]  ra  b*|  rx* 

.c  dj  Ly- 


(1.3) 


£xa  +  yc  xb  +  ydj  £xj 

x  (xa  +  yc)  +  y(xb  +  yd)j 
x  2  a  +  xy(c  +  b)  +  y*dj 


which  means  that  the  system  consisting  of  the  one  equation 

xla  +  xy(c  +  b)  +  ytd  =  e 

is  represented  by  the  matrix  equation  (1.3).  These  two 
examples  show  that  matrix  equations  warrant  attention  if  for 
no  other  reason  than  they  can  help  solve  systems  of 
equations  arising  in  practical  problems.  Indeed,  given  the 
nonlinear  optimization  challenge 


Maximize  /(x)  subject  to 


<1.4) 


equality  constraint x 
equality  constraint 

the  constraint  set  may  be  expressable  as  a  matrix  equation. 

If  so,  this  equation  can  be  solved  to  identify  the  feasible 

region  determined  by  the  constraints  and  thence  the  optima. 

This  approach  is  an  alternative  to  solving  (1.4)  using,  say, 

Lagrangian  multipliers  and  partial  differential  calculus. 

In  examples  (1.2)  and  <1.4),  the  variables  are  often 

taken  to  represent  real  or  complex  numbers.  However,  there 

is  nothing  stopping  the  variables  from  assuming  functional 

values.  For  instance  (1.2)  may  assume  the  form 

a(u,v)x  +  b(u,v)y  ■  c(u,v)  (1.5) 

d(u,v)x  +  e(u,v)y  *  f(u,v) 

where  a,b,c,d,e,f  are  functions  of  the  parameters  u  and  v. 
Though  the  solutions  x  and  y  of  (1.5)  wouldn't  necessarily 
have  numeric  values,  this  is  all  right  because  all  other 
variables  are  functions:  x  and  y  will  likely  take  on 
functional  forms.  To  illustrate. 


[u2v 

u]  [u  +  v] 

fu3v(u  +  v)  +  u(u  -  v)j 

L  uv 

VJ  1m  -  vj 

L  uv(u  +  v)  +  v(u  -  v)J 

■  fu3v  +  u3v3  +  u2  -  uv] 
lu2V  +  uv3  +  uv  -  V3J 

and  so  a  solution  to  the  matrix  equation 


ru3v 

ul  Tx  (u, v)l 

= 

r  U3v  +  U3V=  +  U3  -  uv 1 

Luv 

vJ  Ly<uiv>J 

Lu3V  +  uv3  +  uv  -  V3 J 

is  x(u,v)  *  u  +  v  and  y(u,v)  =  u  -  v.  As  given  in  (1.1), 

this  paper  will  address  matrix  equations  whose  entries 

1.3 


assume  the  -functional  form  of  a  polynomial  in  several 
variables. 

Another  area  where  matrix  equations  appear  (where  the 
matrices'  entries  are  multivariate  functions)  is  control 
theory,  the  science  addressing  the  orbital  natures  of 
objects  in  space.  Matrix  equations  frequently  arise  from 
systems  of  differential  equations  which  describe  a 
satellite's  orbit.  For  example,  given  X(t)  a  columnar 
matrix  whose  entries  (each  describing  an  aspect  of  a 
satellite's  orbit)  are  functions  in  the  parameter  t,  a 
theorem  from  control  theory  states  that  the  equilibrium  of 
the  system  of  differential  equations 

dX 

—  =  FX 

dt 

is  asymptotically  stable  if  the  Lypanov  matrix  equation 

F"rP  +  PF  ■  -C  (1.6) 

has  a  positive-definite  solution  P  for  any  matrix  C  >  0. 
Though  the  terms  in  this  theorem  won't  be  defined  here  (see 
£4:144-1523),  the  important  item  is  that  a  satellite's 
orbital  stability  depends  on  the  solution  of  a  Lypanov 
matrix  equation  (1.6). 

Matrix  equations  can  also  describe  the  path  taken  by  an 
X-ray  passing  through  matter.  This  is  a  concern  of  the 
medical  community  since  the  quality  of  CAT  scan  pictures 
depends  on  the  manner  in  which  the  X-rays  penetrate  the 
skull.  A  good  description  of  the  path  may  allow  finding  a 


brain  tumor,  whereas  a  poorly  chosen  angle  of  entry  may  not 
R.  Vasudevan  C323  develops  higher  order  Riccati  equations 


describing  how  beam  particles  scatter  upon  hitting  matter, 
e.g.  the  skull.  Simply  put,  scattering  is  modelled  by 
including  higher  powers  o-f  a  solution  matrix  X  in  a  Riccati 
equation.  For  instance,  scattering  may  be  roughly  described 
by  the  second  order  Riccati  equation 

XDX  +  AX  XB  +  C  -  0 

and  more  completely  by  the  third  order  Riccati  equation 

XAXBX  +  XCX  +  DX  +  XE  +  F  =  0 
where  the  coe-f-f icient  matrices  A,B,C,D,E,F  reflect  known 
characteristics  of  the  environment  in  which  the  X-rays 
behave.  Higher  powers  of  X  will  more  accurately  describe 
scattering.  Bellman  and  Vasudevan  C23  describe  techniques 
reducing  a  given  Riccati  equation  to  one  of  lower  form. 

This  reduction  ( quasi -1  inear iz at ion)  leads  to  a  series  which 
converges  to  the  actual  solution  of  the  original  equation. 

Literature  Review.  The  literature  supporting  this 
paper  addresses  three  topics: 

I  generalized  inverses  of  matrices 
II  algebraic  structures  of  matrices  over  fields  and 
rings 

III  the  Lypanov  and  Riccati  equations 

Category  I.  An  understanding  of  matrix 
generalized  inverses  launched  the  research  behind  this 
paper.  Matrix  generalized  inverses  encompass  the 
traditional  notion  of  an  inverse  by  assigning  these  to 
non-square  matrices.  The  theory  is  well  developed  for 


matrices  whose  entries  are  complex  numbers*  and  is  described 
below. 

The  physicist  R.  Penrose  C28D  proved  in  1955  that  -for 
any  matrix  A  whose  entries  are  complex  numbers,  there  exists 
matrices  U,Y,Z,W  having  complex  entries  such  that: 

1.  AUA  *  A  (1.7) 

2.  YAY  *  Y 

3.  AZ  equals  its  conjugate  transpose  (i.e. ,  AZ  is 
Her mi tian) 

4.  WA  is  Hermitian 

The  matrices  U,Y,Z,W  are  known  as  the  "generalized  inverses" 
o-f  A.  Penrose  also  showed  the  existence  o-f  a  matrix  X  which 
simultaneously  satisfied  conditions  (1)  thru  (3),  and  a  Y 
simultaneously  satisfying  (1),  (2)  and  (4).  Penrose  also 
showed  that  if  a  matrix  X  satisfied  all  conditions  (given  a 
matrix  A),  then  X  was  unique.  From  the  traditional  matrix 
algebra  viewpoint,  this  unique  X  corresponds  to  the  familiar 
inverse  of  a  square  matrix  A  whose  determinant  is  not  equal 
to  zero.  Though  Penrose's  work  was  original,  he  was  unaware 
that  the  mathematician  E.H.  Moore  [241  had  proved  (using 
abstract  algebra)  the  existence  of  generalized  inverses  for 
arbitrary  rings  more  than  three  decades  earlier. 

Rao  and  Mitra's  textbook  [291  is  devoted  to  the  study 
of  generalized  inverses,  and  is  a  frequently  cited  authority 
in  this  area.  The  first  two  chapters  of  C63  present 
techniques  (with  the  justifying  theorems)  generating  various 
types  of  generalized  inverses  for  a  given  matrix  A,  while 
Captain  Craig  Murray,  AFIT  Class  GCS  85D,  has  recently 


completed  the  programming  o-f  these  techniques.  Two  papers, 
L121  and  C213,  use  generalized  inverses  to  prove  theorems 
about  special  types  o-f  matrix  equations  whose  entries  are 
complex  numbers. 

The  power  o-f  generalized  inverses  weakens  as  one  moves 
-from  matrices  with  complex  entries  to  those  with 
multivariate  entries,  the  subject  o-f  this  paper.  This  is 
because  the  theory  in  the  latter  is  still  very  young. 

Sontag  C31D  discusses  the  existence  o-f  some  generalized 
inverses  -for  special  matrices,  while  Jones  Cl  11  uses  these 
generalized  inverses  to  solve  certain  matrix  equations. 

Category  II.  Most  of  this  paper's  research  dealt 
with  the  decisive  influence  on  the  existence  and  form  of 
solutions  to  matrix  equations:  the  algebraic  structure  of 
matrices.  An  illustration  of  the  importance  of  algebraic 
structure  is  readily  given:  the  solution  of  the  equation 

x  +  7  *  5  (1.8) 

cannot  be  found  among  the  set  of  positive  integers  Z~*  (even 
though  the  coefficients  7  and  5  are  positive  integers) 
because  its  structure  doesn't  include  negative  numbers. 
Since  -2  is  the  solution,  (1.8)  must  be  placed  in  a  larger 
algebraic  structure  in  order  for  a  solution  to  exist.  Thus 
one  moves  from  Z"-  to  the  set  of  all  integers  Z.  Though  Z 
has  a  larger  algebraic  structure  than  1*  it  still  isn't 
large  enough  to  contain  the  solution  of  the  equation 


even  though  the  coefficients  2, -8, 7  are  found  in  Z.  It  now 
becomes  necessary  to  move  to  a  larger  structure  (the  set  of 
rational  numbers)  to  find  the  solution  to  (1.9),  15/2. 
Likewise,  to  solve  the  equations  x3  =  2  and  x3  +  1  =  0  one 
must  move  to  the  larger  algebraic  structures  found  in  the 
irrational  and  complex  numbers  respectively. 

Matrices  too  lend  themselves  to  different  structures. 
However,  not  only  must  the  form  of  a  matrix  be  addressed 
(e.g. ,  whether  it's  invertible,  whether  it's  similar  to 
another  matrix),  but  the  structure  in  which  individual 
entries  are  found  must  also  be  considered.  For  instance, 
decomposing  a  matrix  into  a  product  of  other  matrices  may 
depend  upon  whether  or  not  an  entry  is  irreducible  in  its 
own  setting.  A  trivial  example  is  decomposing  the  1X1 
matrix 

Cx3  +  x  -  13  (1.10) 

into  a  product  of  two  other  1X1  matrices  each  having 
polynomial  entries  whose  coefficients  are  integer,  as  is  the 
entry  in  (1.10).  Such  a  decomposition  is  impossible  because 
x3  +  x  -  1  cannot  be  expressed  as  a  product  of  two 
polynomials  with  integer  coefficients.  Indeed,  the  roots  of 
this  equation  are  (-1  t  *5)  /  2  ,  which  are  irrational. 

In  general,  this  paper  allows  matrix  entries  to  be 
multivariate  polynomials  (having  complex  coefficients,  e.g., 

x3  +  xy,  xy3z  -  <2  +  3i)xayz),  who  in  turn  have  complex 
algebraic  structures.  In  fact,  a  solid  understanding  of 


abstract  algebra  would  have  helped  make  several  articles 


intelligible.  Though  Fraleigh's  text  C73  was  -found  to  be  an 
excellent  primer  on  abstract  algebra,  time  did  not  allow 
sufficient  study  of  a  subject  key  to  this  paper.  Wang  C333 
gives  an  algorithm  for  irreducible  factoring  of  multivariate 
polynomials  having  coefficients  in  an  arbitrary  algebraic 
number  field.  McClellan  C233  presents  methods  for  solving 
systems  of  equations  involving  univariate  polynomials  with 
rational  coefficients  (this  paper  highlights  his  doctoral 
dissertation).  Two  papers  from  the  early  twentieth  century, 
C33  and  E271,  discuss  the  form  of  factors  between 
polynomials. 

The  algebraic  structures  of  matrices  was  addressed  by 
Frost  and  Storey  C81  and  Lee  and  Zak  C203.  These  special 
structures,  called  Smith  Forms,  are  discussed  in  the  next 
chapter.  Unfortunately,  only  the  Smith  Forms  of  matrices 
having  bivariate  polynomial  entries  can  be  addressed:  the 
Smith  Form  involving  multivariate  entries  remains  an 
unresolved  issue  in  mathematics. 

Category  III.  Working  from  a  base  rooted  in 
generalized  inverses  and  recognizing  the  critical  role  of 
algebraic  structures  on  matrices  (and  again  among  their 
entries) ,  a  study  of  matrix  equations  can  proceed.  Though 
much  has  been  published  on  finding  solutions  to  systems  of 
equations,  many  authors  succeed  in  finding  only  particular 
solutions.  After  all,  finding  the  general  forms  of  all 


solutions  may  have  either  been  a  far  too  ponderous  task  or 
the  mathematical  approach  proved  elusive.  In  any  case, 
viewing  a  system  of  multivariate  equations  as  a  single 
matrix  equation  (having  a  solution  in  its  own  right)  is  not 
commonplace.  Depending  on  the  nature  of  the  system,  the 
entire  solution  set  may  be  found  by  solving  a  matrix 
equation  representing  the  system. 

Work  has  been  done  on  matrix  equations  having  complex 
entries.  Roth  £301  (a  frequently  cited  paper)  wrote  on  the 
Lypanov  matrix  equation  in  1952,  and  had  his  work  extended 
in  1972  by  Jones  £121  who  also  addressed  the  second  order 
Riccati  equation,  using  generalized  inverses  to  identify 
solutions.  Morris  and  Odell  £251  attempted  to  find  the 
common  solutions  to  a  set  of  linear  matrix  equations 
<A*x  *  B*}  by  using  generalized  inverses.  Lancaster  £181 
provides  several  approaches  (none  using  generalized 
inverses)  to  salving  J:'B  matrix  equation 

JArXB*  *  C 

As  matrix  equations  assume  multivariate  polynomial 
entries,  the  constraints  cited  in  the  previous  section 
seriously  handicap  the  search  for  solutions.  As  a  result, 
pioneering  work  in  these  equations  is  still  underway. 

The  papers  £113,  £163  and  £173  give  recent  results  connected 
with  the  Lypanov  and  second/third  order  Riccati  equations. 
Collectively,  these  papers  address  the  roles  of  generalized 
inverses  and  matrix  algebraic  structures  in  the  search  for 


II 


Foundations 


This  chapter  lays  the  groundwork  -for  later 

developments.  Because  this  paper  seeks  solutions  to  matrix 

equations,  the  latter's  nature  will  be  addressed  -first. 

Matrix  Equations.  A  matrix  equation  is  an  algebraic 

expression  whose  arguments  involve  matrices.  Matrix 

equations  are  -frequently  used  to  represent  a  kindred  system 

o-f  equations.  For  instance,  the  system  of  equations 

2x  +  3y  =  7 
4x  -  2y  =  1 

is  represented  by  the  matrix  equation 

ii  .;i  p  -  [3 

which  has  the  familiar  form 


where 


AX  ■  b 


(2.1) 


The  reason  for  expressing  a  system  of  equations  in  matrix 
algrebraic  form  is,  of  course,  so  that  the  machinery  of 
matrix  algebra  can  be  used  to  find  a  solution  to  the  given 
initial  system. 

Other  systems  of  equations  have  matrix  representations. 
For  instance,  the  system 


x  +  2u  —  5 
3x  +  4u  =  7 
9x  +  lOu  =  23 
y  +  2v  =  6 
3y  +  4v  =  8 
9y  +  lOv  *  24 
2.  1 


has  the  matrix  -form 
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which  in  turn  has  the  form  (2.1)  where 

A  *  [i  2<]  x  “  [:  V  b 

L  9  10-1 


6 

8 

24 


Thus  in  viewing  <2.1),  it  isn't  necessary  that  b  and  X  be 
one-columnar  matrices  and  A  square,  as  is  traditionally 
presented. 


Lvoanov  Equation.  The  linear  matrix  equation 

(2.1)  leads  to  the  slightly  more  complex  equation 

AX  +  XB  ■  C  (2.2) 

which  is  known  as  the  Lypanov  equation.  In  order  -for  (2.2) 

to  be  de-fined,  the  matrices  A,B  must  be  square  (though 

not  necessarily  of  the  same  size)  and  the  matrices  X,C  must 

share  the  same  size  (though  they  may  be  rectangular).  To 

see  why  this  is  true,  let  the  following  hold: 

MATRIX  A  X  B 

SIZE  a  X  b  cXd  eXf 

For  AX  to  be  defined,  b  must  equal  c:  b  =  c.  For  XB  to  be 

defined,  d  »  e.  AX  is  aXd  and  XB  is  c  X  f.  For 

AX  +  XB  to  be  defined,  a  ■  c  and  d  =  f.  Since  b  =  c  and 

a  ■  c,  a  =  b  and  A  is  square  c  X  c.  Similarly  for  B:  d  =  e 

and  d  -  f  and  so  e  *  f  which  makes  B  square  d  X  d.  It 

follows  that  AX  +  XB  is  c  X  d  which  is  the  size  of  X. 

To  illustrate  (2.2),  the  system 


(2.3) 
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has  the  Lypanov  -form 


(A  X  +  X  B) 
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course. 

(2.3) 

can 

be 

expressed 

as  (2.1) 

by  combining 

terms,  thereby  creating  six  equations  in  six  unknowns. 
However,  it  will  be  seen  (Chapter  IV)  that  by  solving  the 
Lypanov  equation  instead  of  its  equivalent  linear  form 
(2.1),  solutions  of  a  higher  order  matrix  equation  can  be 
derived,  one  in  which  the  Lypanov  is  embedded.  This  higher 
equation  is  known  as  the  second-order  Riccati  equation. 

Riccati  Equation  and  Higher.  In  high  school 
algebra,  a  student  moves  from  solving  the  simple  linear 
equation  bx  =  c  to  the  higher  order  quadratic  equation 
ax2  +  bx  =  c  (in  which  is  embedded  the  linear  equation 
bx  »  c).  In  the  present  setting,  the  former  equation  is 
likened  to  the  Lypanov  matrix  equation  and  the  second-order 
Riccati  matrix  equation  to  the  quadratic.  The  Riccati 
equation  has  the  form 


XDX  +  AX  +  XB  +  C  =  0 


(2.4) 


For  (2.4)  to  be  defined,  A  and  B  must  be  square  (possibly 
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not  the  same  size),  D  must  have  the  size  of  XT,  and  X,C 
share  the  same  size.  The  argument  for  this  follows  that  of 
the  Lypanov  equation. 

To  illustrate  (2.4),  the  Riccati  form  of  the  system 

(2.5) 
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As  with  the  Lypanov  equation,  there  may  exist  many  solutions 
to  the  second-order  Riccati  equation.  In  fact  it  will  be 
seen  in  Chapters  III  and  IV  that  the  solutions  of  (2.4)  will 
satisfy  pairs  of  equations  determined  by  the  matrices 
A,B,C,D  in  (2.4) . 

As  it  is  an  easy  matter  to  create  a  polynomial  of  a 
given  degree  (e.g.,  3x=  +  3x2  -7),  so  is  it  also  to  create 
a  higher  order  matrix  equation.  For  instance,  the  third 
order  Riccati  equation 

XAXBX  +  XDX  +  AX  +  XB  +  C  =  0  (2.6) 

or  the  related  matrix  equation 

EXAXBXF  +  GXDXH  +  AXJ  +  KXB  +  C  =  0 


where  the  matrices  not  equal  to  X  are  known  and  so  can  be 

2.4 


thought  o-f  as  coefficients  of  a  polynomial  in  x.  However, 
the  non-commutativity  of  matrix  multiplication  gives  the 
placement  o-f  each  coefficient  matrix  a  decisive  influence  on 
the  solution. 

A  natural  question  that  arises  is  whether  or  not  a 
given  system  has  an  accompanying  matrix  form.  At  present, 
the  only  known  way  to  answer  this  question  is  to  experiment, 
an  understandably  unattractive  task.  However,  if  a  matrix 
representation  is  derived  (or  stumbled  upon),  the  concepts 
presented  in  this  paper  will  help  to  identify  the  solutions 
to  the  given  matrix  equation,  and  thus  to  its  underlying 
system  of  equations. 

Matrices  Over  a  Ring.  Two  items  determine  the  nature 
of  a  matrix:  its  size  and  entries.  Concerning  the  latter, 
a  matrix  is  said  to  be  "over  a  set  S"  if  and  only  if  its 
entries  come  from  the  set  S.  For  example  if  S  =  Cl,2,3> 
then  the  1X4  matrix 

A  ■  [3112] 

is  "over  S"  because  each  of  its  entries  belong  to  S.  However, 
the  2X3  matrix 


B  =  [  1  2 
L3  2 

is  not  "over  S"  because  B<2,3)  =  0,  which  is  not  a  member  of 
S.  If  however,  S  =  fthe  real  numbers)  ,  then  B  is  "over  S". 
Likewise,  if 


(2.7) 


2.5 


A  -  f2  +  3i  7  ]  (2.8) 

L  9  8  —  7i  J 

then  A  is  not  a  matrix  over  the  real  numbers,  but  is  a 
matrix  over  the  complex  numbers.  In  this  paper,  the  set  of 
complex  numbers  will  be  given  by  the  symbol  <t  (see  the  "List 
of  Symbols"). 

These  examples  lead  to  a  notation  used  throughout  this 
paper:  given  a  set  S,  the  set  of  all  m  X  n  matrices  whose 

entries  are  in  S  is  symbolically  given  by 

S'"*"  (2.9) 

If  A  is  an  m  X  n  matrix  whose  entries  are  in  S,  then  the 
statement  "A  is  an  element  of  SmXri "  is  represented  by 

A  (  S'"*" 

Example:  From  (2.8),  A  6  $zxa!  because  A  is  a  2  X  2 

matrix  whose  entries  are  complex.  However  from  (2.7),  B  is 
not  in  <t3*3  because  Bisa2X3,  despite  the  fact  that  the 
entries  of  B  are  elements  of  the  complex  numbers  (the  set  of 
real  numbers  is  a  subset  of  the  complex  numbers).  However, 

B  6  «=*3  ■ 

The  concept  of  "overness"  will  now  be  extended.  Let 

S(x)  =  {all  polynomials  in  the  variable  x 

with  coefficients  in  the  set  S>  (2.10) 

Thus  if  S  =  {1,2,31,  then 

f(x)  »  xz  +  2x  +  3  f  S(x) 

because  f(x)  is  a  polynomial  in  x  with  coefficients  in  S. 
However 


f (y)  =  y3  +  2y  +  3 


is  not  in  S(x)  because  f (y)  is  a  polynomial  in  y.  Likewise, 

f(x)  ■  7x2  +  2x  +  3 

is  not  in  S(x)  because  the  coefficient  of  x2  is  seven,  which 

s  not  found  in  S  =  <1,2,31.  It  follows  that  <t(x)  is  the 

set  of  polynomials  (in  the  variable  x)  with  complex 

coefficients. 

Similar  to  (2.10),  let 

S(x,y)  =  {polynomials  in  the  variables  x,y 

whose  coefficients  are  in  the  set  S> 

Thus  if  S  =  <4,5,2>,  then 

f(x,y)  =  4x  10ya  +  2x2y  +  5x  +  5y  *  S(x,y> 

because  f(x,y)  is  a  polynomial  in  x,y  with  coefficients  in 

S.  It  follows  that  <t(x,y)  is  the  set  of  polynomials  in  the 

variables  x,y  with  complex  coefficients. 

In  general,  let  r  represent  the  n-tuple  (xi,xa,...,x„), 

Then  given  a  set  S,  S(r)  is  defined  to  be  the  set  of  all 

polynomials  (in  the  variables  found  in  r)  with  coefficients 

in  the  set  S. 

Ex ample:  If  n  =  (x,y,z,v,w>,  then  4(r)  is  the  set  of 

all  polynomials  in  the  variables  x,y,z,v,w  with  complex 
coefficients.  Thus 

f(x,y,z,v,w)  =  (4+5i)x2  +  xyzvw  -  (2+12i)w'7  t(r). 

However  if  S  is  the  set  of  real  numbers,  then  f(x,y,z,v,w) 
i s  not  in  S ( r ) .  ■ 

Now  to  extend  notation  to  matrices.  Define 

S"”<n((’)  =  <m  X  n  matrices  whose  (2.11) 

entries  are  in  S(r)> 


Thus  if  S  ■  fl,2,4,8,9>  and  r  =  (x,y,z)  and 

A  =  r  x2yz  +  9z  x  1 

L  2x  4yz  +  8z2J 

then  A  t  S2x2(p)  because  A  is  2  X  2,  and  each  entry  of  A  is 
a  polynomial  in  p  *  (x,y,z)  with  coefficients  in  S.  However 
if  S  =*  Cl,2,4,8>,  then  A  would  no  longer  be  in  S3x2(p) 
because  a  coefficient  of  A(l,l)  isn't  in  S  (namely,  9). 

From  (2.11),  <tmXr,(r)  is  the  set  of  all  m  X  n  matrices 
whose  entries  are  polynomials  in  p  having  complex 
coefficients. 

With  this  final  definition,  the  object  of  this  paper 
can  be  stated:  find  the  solution  X  fc  Q.f  ^he 

third-order  Riccati  equation 

XAXBX  +  XCX  +  DX  +  XE  +  F  ■  0  (2. 12) 

given  r  the  tuple  of  indetermi nates  in  the  polynomial 
entries  of  A,B,C,D,E,F  (■  $**XJ(p)  for  k,j,m,n  of  appropriate 
size. 

One  may  feel  that  in  (2.12),  mandating  X  k  <t(P)  is 
unduly  restrictive.  After  all,  forcing  the  solutions  of 
x2  -  2  =  0  (a  polynomial  in  x  over  the  integers,  i.e. ,  a 
polynomial  with  integer  coefficients)  to  again  be  a 
polynomial  over  the  integers  isn't  possible,  since  +2  is 
irrational.  Also,  given  the  matrix  equation  AX  3  b  with 
the  entries  of  A,b  integer  and  A  square  and  invertible,  it 
may  not  be  possible  to  have  the  solution  X  ~  A~*b  to  again 
have  integer  entries.  Why  then  require  X  in  (2.12)  to  have 
entries  in  <t(D  when  (like  these  two  analogies)  it  may  be 


necessary  to  leave  <t(n)?  The  argument  is  well  posed,  since 
a  given  matrix  equation  may  indeed  not  have  a  solution  in  a 
given  space.  In  -fact,  practical  solutions  could  be  missed 
by  restricting  solutions  for  theoretical  reasons. 

The  reasons  for  mandating  X  (■  t"*n(p)  are  in  part 
pragmatic.  First  of  all,  polynomial  entries  don't  lend 
themselves  to  singularities:  X  would  be  defined  at  all 
values  of  r.  This  is  important  because  given  an  entry  that 
is  a  rational  function  (e.g.,  1/x),  there  may  exist  a 
singularity  or  discontinuity  at  a  value  which  in  fact  does 
have  physical  significance.  Because  of  the  rational  entry, 
the  form  of  the  solution  may  prove  untenable.  Secondly, 
many  functions  can  be  approximated  by  a  series  of 
polynomials,  e.g.,  Chebysev  polynomials  [5:2391.  In  fact, 
computational  considerations  may  force  the  analyst  to  use 
polynomial  approximations  of  a  function.  Restricting  the 
entries  of  X  to  <t(r)  would  help  generate  these  polynomials 
to  be  used  in  approximating  the  solutions  to  a  phenomena. 

The  driving  force,  however,  for  requiring  X  (  tmXn(r) 
is  theoretical:  is  it  possible  to  find  some  solutions  in 
the  same  space  as  the  coefficient  matrices  in  a  given  matrix 
equation?  If  so,  how  would  the  solutions  be  obtained?  It 
may  very  well  be  that  solutions  thus  found  may  shed  light  on 
the  nature  of  solutions  which  lie  outside  the  required 
space. 

A  word  needs  to  be  said  on  the  setting  for  addition  and 


multiplication  on  matrices  over  4(f).  Though  these 
operations  seem  to  hardly  merit  concern,  there  is 
in  -fact  an  extremely  important  algebraic  structure  behind 
them.  Some  of  the  current  literature  on  matrix  equations 
refers  to  this  "ring"  structure  o-f  4mX"(r>. 

In  the  broad  setting,  a  set  of  objects  S  (e.g. , 
matrices)  is  given,  and  two  operations,  +  and  #,  are  de-fined 
(e.g..,  matrix  addition  and  multiplication)  between  a  pair  o-f 
elements  o-f  S.  These  operations  have  the  property  that  i-f 
a,b  (  S  then:  a+b,a*b£S  and  a  +  b  ,  a  *  b  assume 
unique  values  (e.g.,  a  +  b  has  one  value  and  one  value 
only).  Operations  with  these  closure  and  uniqueness 
properties  are  called  binary  operations.  Depending  on  the 
nature  o-f  the  operations,  an  algebraic  structure  (denoted  by 
<S,  +  ,*>)  is  de-fined  on  the  set  S  taken  together  with  the 
operations.  For  S  the  set  of  matrices  over  <t(r)  and  the 
customary  operations  of  matrix  addition  and  subtraction, 
<S,+,*>  is  known  as  a  ring. 

A  ring  <R,+,*>  is  a  set  R  together  with  two  binary 
operations  +  and  *  defined  on  R  such  that  the  following 
axioms  are  satisfied  C7:1953: 

Rl.  <R,+>  is  an  abelian  group. 

R2.  *  is  associative. 

R3.  For  all  a,b,c  fc  R, 

a*(b+c)  *  (a*b)  +  (a*c) 
and 

(a+b)*c  =  (a*c)  +  (b*c) 

It  is  an  easy  matter  to  confirm  that  the  operations  of 

matrix  addition  and  multiplication  satisfy  the  ring  axioms. 

2.  10 


x..  This  chapter 


Rank  and  Determinant  of  a  liatr: 
closes  by  addressing  two  -familiar  properties  of  a  matrix: 
its  determinant  and  rank.  Both  will  be  cited  throughout 
this  paper. 

The  rank  of  an  m  X  n  matrix  A  (denoted  by  rank (A))  is 
the  number  of  linearly  independent  columns  in  A.  Chapter  III 
will  show  that  rank (A)  affects  the  solution  of  the  matrix 
equation  AX  *  b,  where  A,X,b  are  matrices  over  $(p).  In 
finding  rank (A) ,  the  determinant  of  a  matrix  will  come  into 
play. 

Despite  the  many  techniques  which  find  the  determinant 
of  a  square  matrix,  later  proofs  will  refer  to  the  first 
historical  definition  of  a  determinant,  computationally 
inefficient  though  it  is.  This  definition,  later  referred 
to  as  the  "permutational "  form  of  a  determinant,  derives  its 
name  from  its  indexing  (of  matrix  elements)  on  a  permutation 
of  a  set  C26:B9]. 

Given  the  set  of  integers  S  ■  (1,2,...,N>  a 
permutation  on  S  is  a  one-to-one  and  onto  mapping  (T  from  S 
into  itself,  and  is  represented  by  a  2  X  N  matrix  whose 
top  row  is  S  and  whose  bottom  row  is  a  specific  permutation 
of  S.  For  instance,  the  two  permutations  of  <1,21  are  given 
by 

and  ( 1  2^ 

V  2/ 

For  S  »  <1,2,3,41  one  permutation  is  given  by 


(where  «-(l)  =  2,  T<2>  =  4,  <T<3>  =  3  and  r(4)  =  1)  and 
another  by 

r  -  C  2  3  4) 

\4  3  2  1/ 

where  ffd)  -  4,  f(2)  =  3,  (T<3)  »  2  and  r<4>  =»  1. 

If  A  33  Cah j]  is  an  N  X  N  matrix,  then  the 

determinant  of  A  is  given  by 

det(A)  *  E  Csgn (p>  n j-* (Aj ,p < ^ > ) ]  (2.13) 

pf  S 

where 

S  is  the  set  of  all  permutations  on  C1,2,...N> 

and 

sgn(p)  »  il,  depending  on  the  number  of  transpositions 
which  comprise  the  particular  permutation  p  [7:47-483. 


The  key  to  solving  the  Lypanov  and  Riccati  equations 
involves  -finding  the  solutions  of  the  linear  equation 

AX  -  b  (3.1) 

As  before,  matrices  may  be  non-square  but  entries  must  be  in 

<t  ( P ) . 

Though  the  matrices  of  (3.1)  may  be  rectangular,  one  is 
first  led  to  investigating  properties  of  a  square  matrix 
A  t  $r'Xr,(r) — in  particular,  whether  or  not  A  is  invertible. 
As  in  the  real  case,  if  det(A)  equals  zero,  then  A  has  no 
inverse.  In  general  however,  the  determinant  of  A  is  a 
polynomial  over  r  and  thus  has  an  inverse  except  at  the 
roots  of  det(A).  For  example,  if 

[x  +  y  x 

t  <t3X3(x,y)  (3.2) 

y  *y  J 

then 

det(A)  *  (x  +  y)xy  -  xy  *  xy(x  +  y  -  1) 
and  A  has  an  inverse  except  in  the  root  set  of  det(A) 

C(x,y)  :  x=s0ory*0orx+y=l> 

If  A  must  be  invertible  for  all  values  of  r  (as  will 
soon  be  required) ,  then  det(A)  can  only  be  a  non-zero 
complex  number.  In  this  case  A-1  is  easily  found  using  the 
well  known  formula  126:961 

A-1  =  |A|-1  Cmatri x-of-cof actors  of  ADT 
That  A-1  4  <"Xr,(p)  can  be  seen  by  recalling  that  the 
cofactor  of  an  entry  a4 j  of  A  is  the  determinant  of  a 


submatrix  of  A  (obtained  by  deleting  the  i 'th  row  and  j 'th 
column  of  A).  Because  the  entries  of  A  are  in  $(r),  the 
cofactors  are  also  in  $(r),  and  thus  the  transpose  of  the 
matrix  of  cofactors  is  in  <tr,Xr'(r).  Since  det(A)  is  assumed 
to  be  a  non-zero  complex  number,  so  too  will  !A!~X  which 
implies  A-1  f  4r,xr,(r>.  With  this  in  mind,  a  specific  form 
of  (3.1)  can  be  salved  for  Xs  if  A  6  4"x"(r)  and 
det(A)  6  <  -  f0>  then  there  exists  an  A-1  &  <tr,Xr*(r)  and 
X  =  A~*b  (the  general  setting  for  the  familiar  result 
involving  real  matrices).  Finally,  since  the  entries  of  A-1 
and  b  are  in  <MD,  so  too  are  the  entries  in  their  product 
X  **  A-tb  as  mandated  in  (3.1). 

The  picture  becomes  more  complicated  when  A  doesn't 
have  an  inverse.  This  may  be  because  det(A)  *  0, 
det(A)  t  «(r)  -  t  ,  or  det(A)  isn't  defined,  i.e.  A 
rectangular.  Nevertheless,  (3.1)  may  still  have  a  solution 
and  in  fact,  may  have  infinitely  many.  However,  before  the 
general  solution  of  (3.1)  can  be  derived,  a  more  universal 
setting  for  'inverses'  of  a  matrix  will  be  addressed. 

The  basic  concept  here  involves  elementary  row  and 
column  operations  on  a  matrix  A  f  ♦’BXr,(rl).  Similar  to 
their  counterparts  for  real  matrices,  elementary  row 
(column)  operations  on  A  are  limited  to: 

1.  interchanging  two  rows  (columns) 

2.  adding  a  polynomial  multiple  of  a  row  (column)  to 
another  row  (column) 

3.  multiplying  a  row  (column)  by  a  complex  scalar  not 
equal  to  zero. 


3.2 


As  with  real  matrices,  an  elementary  row  (column)  operation 
on  At  t"*n ( p )  is  equivalent  to  left  (right) 
multiplication  of  A  by  the  corresponding  elementary  row 
(column)  matrix.  For  example,  interchanging  the  rows  of 
(3.2)  can  be  accomplished  by  multiplying  A  on  the  left  by 

[?  i) 

obtained  by  interchanging  the  rows  of  I  2x2.  Another 
example:  multiplying  the  first  column  of  (3.2)  by  x2  and 

adding  it  to  the  second  column  can  be  accomplished  by 
multiplying  A  on  the  right  by 


obtained  by  multiplying  the  first  column  of  I2X3  by  x2  and 
adding  it  to  the  second  column.  Again,  as  with  real 
matrices,  the  determinant  of  an  elementary  matrix  over  $(r) 
is  a  nonzero  complex  number,  and  thus  the  elementary  matrix 
always  has  an  inverse.  It  is  for  this  reason  that  the  third 
elementary  operation  is  limited  to  multiplication  by  a 
scalar  for  if  polynomial  multiples  were  allowed,  the 
determinant  of  the  resulting  elementary  matrix  may  be  a 
polynomial  with  degree  1  1.  The  matrix  would  then  fail  to 
have  an  inverse  at  the  roots  of  its  determinant.  A  square 
matrix  whose  determinant  equals  ±1  is  called  uni modular,  and 
so  is  invertible. 

Both  elementary  row  and  column  operations  are  used  on  A 
to  find  the  general  solution  to  (3.1).  Though  two  different 


cases  must  be  considered — A  does  and  doesn't  have  constant 


rank — both  follow  the  same  approach:  find  uni  modular 
matrices  P  and  Q  that  reduce  A  to  a  form  yielding  the 
general  solution  of  (3.1).  The  two  cases  will  be  analyzed 
separately. 

AX  »  b.  rank (A)  constant.  The  aim  here  is  to  transform 
A  t  $mXr*,  using  elementary  operations,  into  the  m  X  n  matrix 


where  p  is  the  rank  of  A.  The  main  question  thus 
becomes, "What  are  the  unimodular  matrices  P  and  Q  such  that 
PAQ  ■  At?"  Once  P  and  Q  are  found,  the  general  solution  of 
(3.1)  will  follow,  as  this  chapter  will  show.  That  such  P 
and  Q  do  exist  wa  •>  recently  proven  by  E.  Sontag  C313: 

THEOREM  3.0  (Sontaa):  The  following  statements  are 
equivalent  for  a  matrix  A  =  A(r)  over  R  *  $(r): 

a.  There  exists  a  matrix  B  over  R  such  that 
ABA  ■  A  and  BAB  *  B.  B  is  called  the  weak 
generalized  inverse  (or  the  fl ,2>-general ized 
inverse)  of  A. 

b.  There  exists  square  unimodular  matrices  P  and  Q 
over  R  such  that  Ai  =  PAQ. 

c.  As  a  function  of  the  the  complex  variables 
r1  *  (ri,...,rn),  rank  (A)  is  constant. 

Proofs  See  C31D.  ■ 

A  method  for  determining  P  and  Q  has  been  developed  by 
Dr.  John  Jones  Jr.,  AFIT,  and  works  by  keeping  a 
'cumulative'  record  of  all  the  elementary  row  and  column 
operations  performed  in  reducing  A  to  A,.  The  "Jones  ST 


I 


I 


I 

i 


i  *>* 


\ 


> 


> 


0 


Method"  begins  by  -forming  the  matrix 

[AmXr»  |  I  mXm  1 

- -  <3. 4) 

InXn  1  OnXm  " 

Subsequent  elementary  operations  are  done  on  A3  until  A 
becomes  A*.  Once  this  point  is  reached,  the  matrix 
occupying  the  ImXm  block  is  P,  and  the  matrix  occupying  the 
Inxr>  block  is  Q.  That  these  two  blocks  do  indeed  contain  P 
and  Q  is  easily  seen:  P  and  Q  will  be  the  respective 
products  of  the  elementary  row  and  column  transformations 
done  in  reducing  A  to  Ax;  since  ImXm  is  row  adjacent  to  A,  P 
will  be  the  left  multiplier  of  A,  and  since  Inxn  is  column 
adjacent  to  A,  Q  will  be  A's  right  multiplier.  Also,  both  P 
and  Q  are  invertible,  since  their  determinants  are  nonzero 
complex.  To  see  this,  consider  P  ■  n  R*  where  the  's 
are  the  elementary  row  transf ormations  done  in  reducing  A  to 
Ai,  It  follows  that 


det  <P)  -  det ( n  R* )  >  n  det (R* ) 

Since  each  det <R» )  is  nonzero  complex  (the  nature  of  an 
elementary  matrix),  so  will  their  product  be,  and  thus 
det(P).  A  similar  approach  holds  for  Q. 

Example:  Reduce  to  the  form  of  (3.3)  the  matrix 

'  1  x  y 

A  - 

.0  0 


I] 


t  t2X3(x,y> 


Begin  by  augmenting  A  by  the  identity  matrices  I3X2  and 
1 3 x 3  to  form  per  (3.4) 


3.5 


Since  det(P) 


1  and  det(Q) 


— 1,  P  and  Q  are  unimodular  and 


hence  invertible.  ■ 

The  above  example  illustrates  a  theorem  helping  to 
identify  the  general  solution  of  (3.1)  for  A  of  constant 
ranks 


THEOREM  3.1;  Let  A  (  tmXn(r)  and  rank (A)  =  p  a 
constant.  If  P  and  Q  are  unimodular  and 


PAQ 


I  pi  X  p  Op  X  n — p 

• 0m — pXp  0 m — p x n — p - 


where 


P  =  [TJ  f  ,  Q  ■  CS  N3  (  tnXn  ( D  , 

T  i  <t,3Xm(r)  ,  M  i  ?  g  f. 

N  6 

then 


rA  Imxml  is  similar  to  f  Ip*p  0p*n-P  ToX™  i 


x  y  x3  +  l 

0  0  0 

1  0 

-x  y  1 

10  0 

0  1  -1 

03X2 

0  0  1 

Now  add  -x  times  column  1  to  column  3  to  obtain: 


X 

y 

1 

1  0 

0 

0 

0 

-xy  1 

1 

0 

-x 

0 

1 

-1 

0ax2 

0 

0 

1 

, 

Switching  columns  1 

and  3: 

1 

y 

X 

1  0 

0 

0 

0 

-xy  1 

-x 

0 

1 

-1 

1 

0 

03X2 

1 

0 

0 

Adding  -y  times  column  1  to 

column  2,  followed 

by  -x  times 

column  1  to 

column 

3  yields: 

i 

» 

1 

0 

0 

Cl  0  3 

»  T 

0 

0 

0 

C-xy  1  3 

>>  M 

-x 

xy 

1  +  X3 

-1 

l  ♦  y 

X 

03X2 

1 

-y 

-x 

A 

A 

A 

A 

\ 

/ 

S 

N 

Since  A  has 

been  reduced  to 

the  form  (3.3),  P 

and  Q  have 

been  -found. 

and  thus  S,T,M,N  of  Theorem  3.1. 

It's  easy  to 

verify  the  uni modularity  of 

P  and  Q.  ■ 

The  groundwork 

is  now  in  place  to  identify  the  general 

form  of  the 

solution  to  (3.1).  The  following 

theorem  is  due 

to  Jones  £111 


THEOREM  5.2;  Let  A  f  Qf  constant  rank  p  and 


b  fc  4mXv<r>.  Then  AX  *  b  has  a  solution  X  t  <tr,Xv(p)  if 
and  only  if  Mb  *  On,-pXv,  M  given  in  Theorem  3.1.  The 
general  form  of  X  is  given  by 

X  -  STb  +  NZ 


S,T,N  as  given  in  Theorem  3.1,  and  Z  f  arbitrary 

PROOFS  Let  P  fc  «mXm(r)  and  Q  €  <t"x"(P)  be 
unimpdular  matrices  such  that 


PAQ  = 


I  pXp 

i-0„ 


QpXn — p 
.  P  Om — p  .  r» — p  . 


as  given  in  Theorem  3.0.  Then  AX  -  b  has  a  solution  X 

iff 

PAX  =  Pb  has  a  solution  X  <by  virture  of  the  unimodularity 
of  P  and  hence  the  existence  of  P_l) 

iff 

PA(QQ-1)X  «  Pb  =  (PAQ) (Q-XX)  has  a  solution  X  (since  Q  is 
uni modular,  there  exists  Q~x> 


(PAQ) Y  = 

Pb  has  a  solution 

iff 

Y  *  Q~XX 

Let 

Y  * 

fWpXv  1  9  P  = 

\JrXm  1  ’ 

Q  =  CSnXp 

l Zn-pXvJ 

p  Xm  J 

(3. 


Nn  X  n  — p  3 


[Tbl 

LMbJ 


(3. 


Then  substituting  for  PAQ,  P  and  Y  in  (3.5): 

ripxp  0pXn-p  i  po  =  m  b 

L  0m— p>.p>  0m— p  ,n~p  j  LzJ  Lmj 

rw  i 

. Om-pXv J 

Since  Y  ■  Q“XX  from  (3.5), 

X  -  QY  -  C5  N3  £Wj 

Note  that  (3.6)  allows  Z  to  be  arbitrary  since  Z  is  zeroed 
out  by  the  second  row  of  PAQ.  (3.6)  also  says  that  if  a 


SW  +  NZ 


STb  +  NZ 


solution  exists  to  AX  =  b  ,  then  Mb  must  equal  0m.pI>v. 
This  mandate  is  referred  to  as  the  'consistency  condition' 


EXAMPLE  (THEOREM  3.2); 

Given 

x  y  x=  +  y  +  1 

x*y  xy3  x3y  +•  xy3  +  xy 

b  -  r  x  +  xy  +  y  +  x3y  +  y=  -J 

Lx=y  +  xy3  +  x  3y3  +  x3y=  +  Xy3 J 

solve  AX  =  b. 

The  example  far  Theorem  3. 1  (which  reduced  A)  found 

that 

T  =  Cl  01  ,  M  =  C-xy  11  ,  S  =  C-x  -1  ID* 
and 


N 


xy  1  +  x3 
1  +  y  x 
-y  -x 


It  is  a  simple  matter  to  show  Mb  =  0  thereby  confirming 
the  existence  of  a  solution  for  AX  =  b.  From  Theorem  3.2 
it  follows  that  the  general  solution  is 

X  =  STb  +  NZ  = 

■  -x3  -  x3y  -  xy  -  x3y  -  xy3]  f  xy  1 

-x  -  xy  -  y  -  x=y  -  y3  I  +  a  1  +  y  + 

.  x  +  xy  +  y  +  x3y  +  y3 J  L  -y  J 

where  a, ft  k  <t(x,y)  ■ 

A  benefit  of  Theorem  3.2  is  that  it  easily  identifies 

the  basis  for  the  kernal  of  the  linear  transf ormation 

represented  by  the  m  X  n  matrix  A,  a  task  generally  quf. te 

tedious.  Recall  that  the  kernal  of  a  linear  transf ormation 

3.  1 1 


represented  by  the  matrix  A  is  the  set 


K«  »  {  X  :  AX  *  OmXl  } 

To  -find  the  basis  of  KA,  one  must  solve  the  equation 

AX  -  Omxi 

which  is  a  special  case  o-f  AX  ■  b  where  b  =*  0mXi  From 
Theorem  3.2,  it  follows  that 

X  *  STb  +  NZ  *  STOmxi  +  NZ  *  NZ 
Thus  the  span  of  the  columns  of  N  is  K«.  However,  since 
Q  =  CS  N3  and  Q  is  uni modular,  the  columns  of  Q  are 
linearly  independent,  and  thus  the  columns  of  N.  Since  the 
columns  of  N  is  a  linearly  independent  spanning  set  of 
it  follows  that  the  columns  of  N  are  the  basis  for  the 
kernal  of  the  linear  transf ormation  represented  by  the 
matrix  A. 

AX  »  b.  rank < A)  not  constant.  Consider  the  following 
matrix  A  f  <t3X3  (x  ,y  ,  z )  which  doesn't  have  constant  rank: 


Rank (A)  i  <0, 1,2,3}  ,  depending  on  the  values  assumed  by 
(x,y,z):  if  x  a  y  a  z  =  0,  then  rank(A)  =  0;  if  x,y,z  =t=  0, 

then  rank (A)  =  3;  if  two  are  0  while  the  third  isn't, 
rank (A)  »  1;  if  two  aren't  zero  while  the  third  is,  then 
rank  (A)  =*  2.  This  section  addresses  the  solution  of  the 
equation  AX  =  b  given  rank(A>  not  constant. 

A  result  similar  to  Theorem  3.2  holds  for  such  a 


matrix  A.  The  approach  is  practically  identical  to  the 


constant  rank  case:  find  uni  modular  matrices  P  and  Q  which 
reduce  A  to  the  form 

[Ipxp  Opxn— p  1 

I  ,  p  >  1  (3.8) 

— p  Xp  Am — pXn  —  pj 

where  p  is  as  large  as  possible  and  ft  fc  jm-pxn-p^^ 
before,  P  and  Q  are  found  using  the  Jones  ST  Method.  If 
such  P  and  Q  are  found  (they  may  not  exist),  the  following 
theorem  applies: 

THEOREM  3.3:  Let  A  A  $mXr,(p>  jf  p  and  Q  are 
uni modular  and 

[Ip»X,3  Opxn-p 

Om-pXp  An-pXn-p 

where 

P  -  JTj  *  tmxm(r)  ,  Q  ■  CS  N1  f  <t"*"(r), 

T  i  <tpXm  (  n )  ,  M  £  f  g  £  <{nXp(P)  ( 

N  fc  ^nXn—  p>  (  |-1 )  t  A  £  Jm-pXP-p(p) 

then 


r; 

I  mXm  *1 

is  similar  to 

^  P  Xp 

^pXn — p 

T  p  Xm  1 

i.  InXn 

OnXm  J 

^m—p X  p 

Am — pXn — |p 

Mr* — p  Xm  1 

•  SnXp 

Nn  X  n — fp 

On  xm  J 

PROOF:  identical  to  that  of  Theorem  3.1.  ■ 

Once  P  and  Q  have  been  found,  the  following  theorem  is 
appl ied: 

THEOREM  3.4:  Let  A  (  <tmXr*(r)  and  b  t  <tmX''(r').  Ther 
AX  «  b  has  a  solution  X  £  4r,x',(r)  if  and  only  if  there 
exists  a  Z  i  <tr'-,3*v(n)  such  that  Mb  -  AZ  M, A  given  in 


Theorem  3.3.  In  this  case,  the  general  form  of  X  is  given 


X 


STb  +  NZ 


■for  all  Z  such  that  Mb  *  AZ  and  S,T,N  given  in  Theorem  3.3. 

PROOF :  Let  P  4  and  Q  4  be 

uni modular  matrices  such  that 


PAG 


•I  =  x.  °-‘ -  I 

•  0(B-p  ,p  ilm-p  fn— p  J 


AX  =  b  has  a  solution  X 


iff 


PAX  =  Pb  has  a  solution  X  (by  virture  of  the  unimodularity 
o-f  P  and  hence  the  existence  of  A-1) 

iff 

PA(QQ-1)X  =  Pb  =  (PAQ) (Q-XX)  has  a  solution  X  (since  Q  is 
uni  modular,  there  exists  Q~x) 


iff 

(PAQ)  Y  =  Pb 

has  a 

solution 

Y  *  Q-x  X 

(3.9) 

Let 

Y  =* 

rw„xv  I , 

p  ~  [t..;.  -I 

l ^n-pXvJ 

L  nm-pxm J 

and 


Q  =  C  S0Xp  NnXn-p 1 

Then  substituting  for  PAQ,  P  and  Y  in  (3.9): 


b 


Since  Y  =  Q~XX, 

X  =  QY  ■  CS  ND  JWJ  -  SW  +  NZ  *  STb  +  NZ 
subject  to 

AZ  =  Mb 


(3. 10) 


(3.11) 


per  (3.10).  This  latter  mandate  is  the  consistency 
condition.  That  Z  <  <tr,-t,Xv(r)  proceeds  from  (3.9):  since 


Q  has  polynomial  entries  and  is  uni  modular,  so  too  is  Q-1. 
Since  X  is  also  required  to  be  in  tr,Xv(r,>  the  product  of 
Q-1  and  X  must  be  in  $"*'•' (r).  Since  Z  is  a  submatrix  of 
Y  -  Q"1X  ,  Z  {  4n-pXv(^  .  ■ 

In  general,  there  is  no  assurance  that  the  consistency 
condition  (3.11)  can  be  met  -for  a  given  equation  AX  =  bs 
though  A  may  be  reduced  to  the  form  (3.8),  there  may  not 
exist  a  Z  OVER  $(r)  satisfying  the  consistency  condition 
AZ  =  Mb.  It  may  be  necessary  to  leave  $(r)  in  order  to  find 
a  Z  satisfying  the  consistency  condition,  an  act  contrary  to 
this  paper's  aim. 

What  then  can  be  said  about  the  solvability  over  4(p) 

of  AX  *  b  when  rank (A)  isn't  constant?  It's  possible  that 

A  isn't  reducible  to  the  form  (3.8).  For  instance  the 

matrix  in  (3.7) s  reducing  A  to  the  form  (3.8)  implies  that 

% 

A  has  a  rank  no  less  than  p  (the  first  p  columns  of  (3.8) 
are  linearly  independent  despite  the  form  of  A).  However, 
the  rank  of  A  may  be  zero  if  x  =  y  =  z  =0  ,  and  so  A  isn't 
reducible  to  the  form  (3.8). 

This  example,  however,  is  symptomatic  of  a  larger 
problem.  Assume  that  a  given  matrix  A  has  been  reduced  as 
much  as  possible  (no  known  method  exists  that  insures 
maximum  reducibility  of  a  matrix  over  $(r>,  r  arbitrary), 
i.e.  the  maximum  value  of  p  has  been  found.  One  is  then 
left  with  finding  a  Z  satisfying  the  consistency  condition. 

A  in  turn  cannot  be  reduced,  since  this  would  violate  the 


3.  15 


assumption  that  A  has  been  reduced  as  far  as  is  possible. 

If  A,  is  square,  det(A)  cannot  be  a  nonzero  complex  number 
(if  it  were,  rank (A)  =  m  -  p  and  A  could  be  reduced  to 
Im-oxm-p,  thus  giving  A  constant  rank,  a  contradiction).  If 
det  (A)  =>=  0  repeated  use  of  Cramer's  rule  may  identify  some 
solutions  Z,  but  these  would  be  local  since  the  roots  of 
det<A>  must  be  omitted.  In  general,  however,  A  will  be 
rectangular . 

Thus  finding  the  general  solution  of  (3.1)  given 
rank (A)  not  constant  ultimately  requires  solving  (3.11) 
where  A  cannot  be  transformed  into  the  form  (3.8).  The 
search  for  Z  may  be  aided  by  matrices  related  to  A,  the 
topic  of  the  next  section. 


SMITH  FORM  OF  A  MATRIX.  A  matrix  A  fc  <t*"Xr,<r)  is  said 
to  be  in  Smith  Form  (SF)  if  and  only  if 


f ! 

• 

0 

0 

• 

f. 

0  j  xn  —  J 

Om- 

J  X  J 

Om— J Xn — J 

(3. 12) 


where  all  off  diagonal  entries  equal  zero  and  f*  is  a  factor 
of  f„-*i  .  For  example 


Tx  0  0  -j 

A  =*  0  x2y  0  i 

LO  0  x  2y3J 


<t3X3(x,y) 


is  in  SF  (3.12)  because  all  off  diagonal  entries  are  zero,  x 
is  a  factor  of  x2y,  and  x2y  is  a  factor  of  x2y=.  Another 


A 


i  <t3XZ  (x  ,y ) 


r  x  0 

0  x  (x  +  1) 

Lo  0 

which  is  in  SF  (3.12)  because  all  off  diagonal  entries  are 
zero,  and  x  is  a  -factor  of  x(x  +  l).  The  solving  of  matrix 
equation  (3.1)  may  be  aided  if  A  can  be  transformed  (using 
unimodular  matrices  P'  and  Q')  into  a  matrix  having  SF.  In 
fact.  Smith  Forms  have  already  appeared:  (3.3)  is  a  special 
case  where  fk-  1  ,  k  <  p  The  main  aim  here  is  to  use 
SFs  to  solve  the  equation  (3.11)  til  =*  Mb  of  the  last 
section. 

To  illustrate,  assume  that  in  solving  JIZ  =  Mb  (derived 
from  some  hypothetical  original  equation  AX  »  b) , 
unimodular  matrices  P'  and  Q'  have  been  found  that  give 
P'AQ'  a  SF.  Thus 

AZ  =  Mb 
P'AZ  =  P  'Mb 
P ' A(Q 'Q ' -l ) Z  =  P 'Mb 
(P'AQ) <Q'-*Z)  =  P  'Mb 

Letting 

Y  =  Q ' _1 Z 
H  =  P 'Mb 

F  *  P'AQ'  (3. 13) 

the  last  equality  becomes 

FY  -  H 

Since  F  is  in  SF,  there  exists  at  most  one  nonzero  entry  in 
each  row  and  coiumn  of  F.  Because  F  and  H  are  known,  some 
entries  in  Y  can  be  found  by  dividing  the  respective  F  and  H 
factors  (if  however  a  given  entry  of  Y  isn't  compatible  with 


two  or  more  equations,  or  if  the  entry  doesn't  turn  out  to 


be  a  polynomial,  then  AZ  =  lib  has  no  solution  Z.  > 


Depending  on  the  structure  of  F,  there  may  be  entries  of  Y 
which  can  assume  arbitrary  values.  Y  is  then  substituted 
into  (3.13)  to  yield  Z  =  Q'Y  ,  thus  solving  (3.11)  and  in 
turn  <3. 1 ) . 

Given  A  {  f  if  there  exists  uni  modular 

matrices  P*  t  <tmXm(r)  and  Q'  t  <t"Xr>(r>)  such  that 
A'  =  P'AQ'  is  in  SF,  then  A  is  said  to  be  equivalent  to  its 
Smith  Form.  It  has  been  proven  11:1883  that  for  every 
A  t  tmXn (x ) ,  A  is  equivalent  to  a  unique  Smith  Form  (an 
algorithm  exists  Cl: 1923  which  finds  this).  However,  such 
is  not  always  the  case  for  multi -dimensional  r.  In  fact, 
Lee  and  Zak  C203  prove  that  a  matrix  A  t  4mXr’(x,y)  is 
equivalent  to  its  SF  iff  a  certain  system  of  linear 
polynomial  equations  has  a  solution.  The  conditions  under 
which  a  matrix  in  three  or  more  variables  is  equivalent  to 
its  Smith  Form  remains  an  open  question. 

Although  the  existence  of  a  SF  for  a  given  matrix  A 
must  be  known  before  a  search  for  it  can  begin,  a  method  is 
needed  which  identifies  the  uni  modular  P'  and  Q'  which 
yield  a  particular  SF  of  A.  As  before,  the  Jones  ST  Method 
can  be  used. 

EXAMPLE:  Find  the  Smith  Form  of  the  matrix  (Frost  and 


Storey  CB3) 


’  s  + 
A  =  0 

0 


0 

s  +  z 
0 

3.  18 


i  <t3X3(s,z> 


^  -*-V 


Next,  add  -<z  +  s)  times  the  third  column  to  the  first 
column  to  get: 


Finally,  interchange  the  first  and  third  columns: 


Since  the  upper  left  hand  matrix  is  in  Smith  Form  (i.e. , 

1  is  a  factor  of  s  +  t  ,  in  turn  a  factor  of  -s<s  +  z)  and 


off  diagonal  entries  are  zero) 


IV  Solutions  to  Lvpanov  and  Riccati  Equations 

This  c.hapter  will  consider  the  solutions  to  the  Lypanov 
matrix  equation 

AX  -  XB  =  C  (4.1) 

and  the  second  order  Riccati  matrix  equation 

XDX  +  AX  +  XB  +  C  =  0  (4.2) 

where  the  entries  of  A,B,C,D  are  in  4(p).  It  will  be  shown 
that  solving  these  two  equations  for  X  hinges  upon  solving' 
equations  of  the  form  AX  —  b  the  subject  of  Chapter  III. 

Lvpanov  Equation.  (4.1)  may  be  solved  using  tensor 
products  or  similarity  transf ormations.  The  former 
approach  will  be  addressed  first  because  it  uses  the  results 
of  Chapter  III  at  the  onset.  The  approach  using  similarity 
transf ormations  will  then  be  presented  as  a  lead-in  to 
solving  the  Riccati  equation  (4.2),  whose  solution  relies 
exclusively  on  these  tranf ormations.  Whichever  method  is 
taken  though,  the  linear  equation  (3.1)  will  demand 
resolution. 

Lvpanov  Equation;  Tensor  Solution.  Multiplication 
between  two  matrices  is  easily  extended  to  that  of  their 
tensor  product.  Given  A  *  Cat j]  the  tensor  product  of  the 
matrices  A  and  B,  represented  by  AaB,  is  given  by  [auB], 
Unlike  matrix  multiplication,  the  number  of  columns  in  A 


doesn't  have  to  equal  the  number  of  rows  in  B 


Ex ample;  If 


ri 

2 

31 

and 

B  =  ri  n 

U 

5 

6  j 

U  2  J 

then 


AaB  =  CaA  jB3  =  flB  2B  3B‘ 

I.4B  5B  6B. 


r  i  i 

2  2 

3  3-1 

1  2 

2  4 

3  6 

4  4 

5  5 

6  6 

L  4  8 

5  10 

6  12  J 

It  follows  that  ifAismXn  and  B  is  r  X  s,  then  AaB  is 
mr  X  ns.  ■ 

The  next  step  in  solving  (4.1)  is  to  identify  the 
solution  of  the  equation 


AXB  =  C 


(4.3) 


where  A*  denotes  the  i 'th  row  of  A  and  A*  its  i 'th 
column  (the  i-jth  entry  of  A  is  given  by  a* j  with  no  comma 
between  subscripts) . 

THEOREM  4.1s  Solving  the  matrix  equation 

AmXnXnXpBpXq  =  CfnXq 

is  equivalent  to  solving  the  equation 

Gu  =  c 


where 

G  «  AaBT  ,  u  =  CXt - X„. «-!-«-  ,  c  =  [Clfl - Cm,r]^ 

PROOF'.  Let  A  =  CaXJ3  ,  B  =  [Bu]  , 

0  ~  Ec  i  j  ]  f  X  —  Ex  i  j  ] 

Since  AXB  =  C  , 


4.2 


w * J  nifrAOj,c 

*  Ax.,-  CXx  .rBj  .C.  .  .  X„.,-Bj  ,e:3'r 

*  Ck [4i i<Ei_  (x j  >  ] 

“  2k  ,i_  <aXkRj_jXki_> 

=  2l_  (suffLjXu.)  +  +  Ei_  (ai„8LJx„u) 

=  aiiEi-(fiLjx1L.)  +  •••  +  a^Ecffi^x^) 

“  c  axiEfixj. .  ...  amCOi  3 

*  C  Xx.. - Xr,.-  ]T 

*  CanBTj,c  ...  alr,BTj 

*  CXi . r . . . Xn ,r 3T  (4.4) 

Since  c  =  CCx . .Cm.,_3T  is  given  and  from  (4.4),  the 
theorem  follows.  ■ 

The  following  theorem  is  due  to  P.  Lancaster  C1B1: 

Theorem  4.2:  Given  the  matrix  equation 

AiXBx  +  A3XB3  =  C  (4.5) 

where  Ax  are  m  X  m  ,  Bx  are  n  X  n  ,  X,C  are  m  X  n  (with  Xr 

and  Cx  the  rows  of  X  and  C  respectively) ,  (4.5)  is  equivalent 

to  solving  the  equation  Gx  =  c  where 

B  ■  AxaBxt  +  A3aB2t 
x  ■  CXx... Xm]T 
c  *  CCx...Cm]T 

Proof :  Similar  to  that  of  Theorem  4.1.  ■ 

The  restriction  on  square  Ax  and  Bx  in  Theorem  4.2  is 
easily  modified  to  include  rectangular  matrices.  With  this 
modification.  Theorem  4.2  allows  for  the  solution  of  the 
Lypanov  equation  (4.1)  which  is  a  special  case  of  (4.5)  with 
Bx  =  I  and  A2  =  -I.  Once  a  given  Lypanov  equation  has 
been  reduced  to  an  equation  of  the  form  AX  -  b  using 
tensor  products,  the  results  of  Chapter  III  can  be  applied  to 
this  latter  equation. 

Example:  Solve  the  Lypanov  equation 


AX  +  XB  ■  C  ■  AX  I  +  IXB 


where 


■  s+z  s  ' 

-  -z  -s’ 

A  = 

,  B  = 

.  z  s-z  . 

.  -z-s  s-z - 

C  = 


s2 

4sz-2z3-3sa 


X  = 


(4.6 


sz 

L2s=-2z=-sz 
‘  Y  z  ' 

.V  w  . 

The  -first  step  is  to  reduce  the  given  Lypanov  equation 
to  the  -form  Gx  *  c  so  that  Theorem  4.2  may  be  used.  Thus 
G  =  A*  IT  +  I *BT  = 

-z  -z— « 


r  s+z  s 
z  s-z 


■1 

o- 

*  1 

O' 

A 

.0 

1. 

+ 

.0 

1- 

* 

_ 

I -S’ 

>-z. 


-s+z 

0 

s  0  ' 

-  — 

z 

-z  -s 

0 

0  -| 

0 

s+z 

0  s 

+ 

— 

5 

5-Z 

0 

0 

z 

0 

s-z  0 

0 

0 

-z 

-z-s 

.  0 

z 

i 

N 

1 

in 

o 

-  0 

0 

-s 

s-z  . 

s 

-z-s 

s 

0 

- 

-s 

2s 

0 

s 

z 

0 

s— 2z 

-z-s 

0  z  -s  -2s-2z  J 
Let  x  =  Cy  z  v  w]'r  and 

c  =  Csz  s2  2s2-2z2-sz  4sz -2z 3-3s2 ] T .  Thus  the 

solution 

■  y  z 

X  ■ 

.  V  w 

to  the  given  Lypanov  equation  is  the  solution  to  the 


linear  equation 


-  s  -z-s  s  0 

-s  2s  0  s 

z  0  s-2z  -z-s 

_  0  z  -s  -2s-2z_ 


rvi 

sz 

Z 

S 

s2 

V 

2s2-2z=-sz 

.  W_ 

_4sz-2z2-3s2 

Depending  on  the  nature  o-f  the  Q  matrix  obtained  using 
tensor  products,  the  given  Lypanov  equation  may  or  may  not 
be  easily  solvable.  For  instance,  if  G  cannot  be  reduced  to 
the  form  (3.3),  one  then  tries  for  form  (3.8).  Here 
however,  finding  the  Z's  satisfying  the  consistency 
condition  (3.11)  flZ  =  Mb  may  be  difficult.  In  any 
event,  the  use  of  tensor  products  creates  a  large  matrix  G, 
thereby  compounding  calculations. 

The  second  approach  to  solving  the  Lypanov  equation 
(4.1)  uses  similarity  transf or mat ions.  Although  this  method 
may  not  explicitly  identify  the  solution  to  (4.1),  it  will 
identify  a  set  of  matrices  within  which  the  solutions  of 
(4.1)  will  be  found. 

Lypanov  Equation;  Solution  via  Similarity 
Transf or mat ions.  The  second  approach  to  solving  (4.1) 
proceeds  from  a  few  simple  observations  unrelated  to  solving 
a  matrix  equation.  Given  an  n  X  n  matrix  X  and  I  =  I„Xn  the 
matrix 


always  has  det(E)  =  1  despite  the  nature  of  X  (this  follows 
from  the  permutational  definition  of  a  determinant  (2.13). 


Since  E  is  uni  modular,  it  has  the  uni  modular  inverse 

4.5 


which  is  verified  below: 


EE'  ■ 


■I  XT  r I  -X i  r 1  +  XO  -IX  +  xn  rl  On 

.0  iJ  LO  IJ  LOI  +  10  -OX  +  i  =  J  Lo  iJ 


(E'E  =  Iznxzn  likewise  -follows). 

The  next  observation  uses  the  matrix 


A  C' 
0  B. 


(4.9) 


where  A,B,C  are  all  nXn  matrices: 


ERE_1  =  [5  5]  [$  g]  [l  -?]  *  ts  CT]  [l  1] 

=  J-A  -AX  +  C  +  XB] 

Lo  B  j 


(4. 10) 


The  above  two  observations  can  now  be  brought  together 
to  solve  (4.1):  the  1-2 'th  entry  of  ERE-1  above  is  itself  a 
Lypanov  equation.  If  it  equalled  zero,  then  X  would  be  a 
solution  to  the  equation 

-AX  +  XB  +  C  *  0  or  AX  -  XB  =  C 
which  is  of  the  form  (4.1).  The  following  theorem  due  to 
Jones  Clll  has  now  been  proven: 

Theorem  4.3:  If  X  is  a  solution  to  the  Lypanov 


equation 


AX  -  XB  =  C 


where  A,B,C,X  f  <t"Xri(r)  then 


,  \  %  *.  *.  % 


R  =  is  similar  to  R'  =  (4.11) 

L  0  BJ  L  0  BJ 

Proof'.  previous  discussion.  ■ 

The  following  more  powerful  theorem  was  proven  by  Roth  C301 

for  matrices  whose  entries  are  complex  numbers: 

Theorem  4.4:  If  A,B,C  t  then  AX  -  XB  ■  C  has 

a  solution  X  £  <tr,x"  if  and  only  if  R  and  R'  in  (4.11)  are 

similar. 

Proofi  See  C301.  ■ 

Although  Theorem  4.3  is  the  core  result  leading  to  a 
solution  of  (4.1),  the  final  assault  upon  (4.1)  will  use 
polynomials  of  a  matrix.  For  this  reason,  the  following 
theorems  are  presented. 

Theorem  4.5:  Let  A,B  be  n  X  n  matrices  such  that  B~l 

exists.  Then  for  n  £  1* 

% 

(BAB-1)"  =  BA^B-1  (4.12) 

Proofi  Inductive  reasoning  will  be  used.  (4.12)  is 
trivially  true  for  n  -  1.  For  0*2, 

(BAB”1)3  =  (BAB~l) (BAB"1)  =  BA (B~lB) AB-1 
=  BA ( I ) AB~ 1  *  BA=B~ 1 

and  (4.12)  is  again  true.  Assume  (4.12)  is  true  for  n  =  k 
Thus 

(BAB"1)1*  a  BA^B-1 

and  so 

(BAB-1)**-1  =  (BAB"1)**  (BAB-1)  =  BA**B-1  (BAB-1 ) 

=  BA**  (B“1B)  AB"1  =  BA**- 1  B“ 1 

and  so  (4.12)  holds  true  for  k  +  1  and  therefore  for  all 


positive  integers.  ■ 


Theorem  4.6:  Let  f(x)  be  a  polynomial  in  the 
variable  x.  Let  A,B  be  n  X  n  matrices  such  that  B-1  exists. 
Then 


f  <BAB~M  =  B  -f  (A)  B-1 

Proofs  Let  f(x)  =  aQ  +  Ek-i.h  aKXK  Then 

f(BAB~M  =  a0I  +  Ea«  (BAB-1  )K  ■  (-from  Theorem  4.5) 
=  a0BB-1  +  Ea^BA^B-1 
*  Ba0B-1  +  B  <  EaKAK)  B~l 
=  B  (a0I  +  EaKAK)  B“l 
■  B  f (A)  B-x  ■ 

Theorem  4.7s  Let 


R 


A,B,C  nXn  matrices.  Then  -for  K  £  Z" 

AK  0  1 


R*<  ■ 


C* 


(4. 13) 


*  an  expression  in  A,B,C. 

Proofs  Inductive  reasoning  will  be  used.  (4.13)  is 
trivially  true  for  i  *  1 ,  where  *  *  B  *  OA  +  IB  +0C 
For  i  =  2  , 


A=  0 

BA+CB  C2 

and  (4.13)  holds  true  again.  Let  (f  13)  be  true  for  n  *  k 
Then 


’A1*  0  i 

a> 

o 

i 

•A**1  0  ' 

.*  C*. 

- 1 

u 

_ I 

. 

.  #  c* 

and  thus  (4.13)  is  true  for  all  positive  integers.  ■ 


Let 


Theorem  4.8; 

R 


■A  0 
-B  C 


A,B,C  nXn  matrices.  If  f(x)  is  a  polynomial  in  x,  then 

r  f  <A>  0  *] 

f  <R>  -  (4. 14) 

L  *  f(C)J 


*  an  expression  in  A,B,C 

Proofi  Inductive  reasoning  will  be  used.  Let 
■f(x)  =30  +  Ek-i  .r.a*x“  For  n  =  1 ,  f  (R)  =  a0  +  atR  = 


[I 

°1 

+  at  f A 

°1 

= 

r  a0 1  +  atA 

0  1 

lo 

IJ 

LB 

CJ 

l  axB 

®o  I  +  a  i  Cj 

f  (A)  0  I 

*  f  (Oj 


where 

*  »  aiB  =  axB  +  OA  +  OC 

and  so  (4.14)  holds. 

For  n  *  2,  f  <R)  =  a0I  +  axR  +  a2R3 


38  a0  T I  Ol  +  ax  TA  01  +  a3  [A3  0  1  (from  Theorem  4.7) 

Lo  iJ  IB  Cj  L*  C3J 

=  J-  a0I  +  axA  +  a2A2  0  ■, 

L  *  aQI  +  aiC  +  a2C2J 

=  ff(A)  O  1 

L  *  f(C)J 

and  (4. 14)  remains  true.  Let  (4. 14)  be  true  for  n  =  v 
If  fi(x)  =  C  a0  +  2k-i,vakxk  ]  and 
f(x)  =  f  i  <x )  +  av+1xv'"1  then 

f(R)  *  [ft (A)  0  1  +  av+lRv+ »■ 

L  *  fi<C)J 

=*  f  ft(A)  0  ]  +  av*i  [A~* 

l  *  fi(C)J  [* 

4.9 


-  -*■  — r-  -A-.  -A  . 


and  so  (4.14)  is  true  for  all  positive  integers.  ■ 


One  more  theorem  remains  to  be  cited  before  the  assault 
on  (4.1)  can  begin.  The  Hami 1 ton-Cay ley  Theorem  will  permit 
the  identification  of  a  pair  of  linear  equations  whose 
common  solutions  will  contain  those  of  (4.1). 

Theorem  4.9  (Hami lton-Cavlev? :  Given  A  an  nXn  matrix 
and  J(y)  =  det (A  -  plnxn)  the  characteristic  equation  of  A 
then 

5(A)  =  0n  Xn 

Proof i  see  Nering  C26: 1003  ■ 

The  tools  are  now  in  place  with  which  to  solve  (4.1). 
The  following  theorem  is  due  to  Jones  Ell  3: 

Theorem  4.10s  Given  the  Lypanov  equation 

AX  -  XB  =  C  AtB,C  fc  4"*" (P)  (4.15) 

let  R  and  R'  be  defined  as  in  (4.11)  and  fA  the 
characteristic  equation  of  A.  If 

r  U  M- 

f« (R)  =  (4.16) 

L  V  N. 

then 

U  +  XV  *  0  (4. 17) 

M  +  XN  -  0 

Thus  a  solution  to  (4.15)  will  be  found  among  the  common 
solutions  of  the  pair  of  equations  (4.17). 


Proofs 


From  the  discussion  preceeding  Theorem  4.3,  it 


was  seen  that 
ERE"1  = 


Vo  i]  [5  g]  [l  1]  -  [5  2] 


=  R  ‘ 


for  X  satisfying  (4.15).  By  substitution, 

f*(ERE”M  =  f^(R') 


(4.  18) 


From  Theorem  4.8, 
f«(R)  = 


ff«(A)  0  1 

1  *  f«(B)J 


’  [l 

From  Theorem  4.6 


0 

f«(B> 


] 


(Hami 1 ton-Cayley  Theorem  4.9) 


f«(ERE~M  »  E  f«(R) 

ii 

H 

1 

Ul 

[l 

XI 

i  ru 

Ml 

ri 

-XT 

II 

C 
<  f 

X 

< 

M+XNl 

N  J 

Lo 

U 

ij 

-?] 

1  L* 

Nj 

Lo 

IJ 

[U+XV 

L  v 


-(U+XV)X  +  (M+XN) 1 
-VX+N 


(4. 19) 


J 


'0 

* 


0 

f*(B), 


(from  (4.18)) 


(4.20) 


Since  two  matrices  are  equal  iff  respective  entries  are 
equal,  the  1-lth  entries  of  (4.19)  and  (4.20)  imply 

U  +  XV  =  0  (4.21) 


and  the  l-2th  entries  imply 

0  =  -(U  +  XV) X  +  (M  +  XN) 

*  OX  +  (M  +  XN)  (from  (4.21)> 

=■  M  +  XN 

Thus  a  solution  of  (4.15)  will  satisfy  the  pair  of  equations 
(4.17).  ■ 

The  example  used  to  illustrate  the  tensor  solution  to 
(4.1)  will  serve  to  illustrate  Theorem  4.10: 


example  Theorem  4.10:  Solve  the  Lypanov  equation 


AX  +  XB  =  C 

*  AX  -  X (-B) 


(4.22) 


where 


z  s-z 


2s2-2z2-sz  4sz-2z2-3s2 


From  (4.11), 


R  =  f  A  C]  -  rs+z 
Lo  -Bj  z 
0 

.  0 


2s2-2z2-sz  4sz-2z2-3s: 


Next,  -find  ,  the  characteristic  equation  o-f  A. 


■f«(p)  *  det  rs  +  z  -  p 


s-z 


*  (s+z-p)  (s-z-JJ)  -  sz 

=  5=  -  sz  -  SP  +  ZS  -Z2  -z  V  -PS  +Pz  +  P2  -  SZ 

=  s2  -  SZ  -  z2  +  p=  -  2s p 

Now  evaluate  -fA(R)  (this  was  done  using  the  program  in 

Appendix  A,  a  BASIC  program  performing  operations  on 
matrices  whose  entries  are  multivariate  polynomials): 


f«(R)  =r0  0  -s2z  +  3s3 

0  0  5sz2  -  5s3  -2z3  +  2s2z 


■ts*  -  2sz 
-sz  -  3s2  +  2z: 


-5s3  +  7s2z  -  2sz2 
2sz2  +  8s3  -  Bs2z 


-3s2  +  2sz 
5s2  -  4sz 


'fU  Ml 
LV  Nj 


According  to  Theorem  4.10,  the  solution  to  (4.22)  will  be 


-found  among  the  common  solutions  to  the  pair  o-f  equations 


U  +  XV  =  0  (4.23) 

M  +  XN  =  0  (4.24) 

Since  U  -  V  =  0,  (4.23)  doesn't  contribute  to  the  solution. 

So  the  search  restricts  itself  to  solutions  of  the  equation 

M  =  -XN  or 


r  -s2z 

+ 

3s3 

-5s3  +  7s=z  -  2sz 

L  5sz2 

— 

5s3  -2z3  +  2s2z 

2sz2  +  8s3  -  8s2: 

=  fw 

x  ■ 

1  r -2s2  +  2sz 

3s2  -  2sz 1 

Ly 

z. 

|  [  sz  +  3s2  -  2z2 

-5s2  +  4szJ 

However,  in  order  to  apply  the  techniques  of 
Chapter  III,  one  must  have  an  equation  of  the  form 

AX  -  b 

It  is  an  easy  matter  to  convert  (4.24)  into  the  form  of 
(3.1)  by  observing  that  M  ■  -XN  implies  MT  =  -NTXT  This 
latter  equation  is  in  the  form  (3.1),  and  the  techniques  of 
Chapter  III  can  be  applied  to  this  transposed  equation.  ■ 

It  is  possible  to  further  restrict  the  set  of  equations 
within  which  the  solution  of  (4.15)  belongs: 

Theorem  4.11:  Given  the  Lypanov  equation 
AX  -  XB  =  C  A , B , C  £  <t"*"(n) 
let  R  and  R'  be  defined  as  in  (4.11)  and  fB  the 
characteri stic  equation  of  B.  If 


fB  (R)  = 


then 


N  -  VX  =  0  (4.25) 

M  -  UX  =  0 


Thus  a  solution  to  (4.15)  will  be  found  among  the  common 


solutions  of  the  pair  of  equations  (4.25). 

Proof :  The  proof  is  similar  to  that  of  Theorem  4.10. 

From  the  discussion  preceeding  Theorem  4.3,  it  was  seen  that 


ERE-1  =  fl  Xl  fA  Cl  fl  -XI  =  fA  01 

Lo  Ij  [o  bJ  Lo  ij  Lo  BJ 

for  X  satisfying  (4.15).  By  substitution, 

f  s (ERE-1 >  =  f * (R ' ) 

From  Theorem  4.8, 


=  R ' 


(4.26) 


f » (R ' )  *  rf»(A) 


0  1 
f B (B) J 


:b(A)  01  (Hami 1 ton-Cay ley  Theorem  4.9)  (4.27) 


From  Theorem  4.6 

f» (ERE-1 )  =  E  f»(R)  E-1  =  fl  X 


U+XV  M+XN 
.  V  N 


fl  xl 

fu 

Ml 

f  I 

-Xl 

Lo  ij 

Lv 

Nj 

lo 

IJ 

l+XN]  fl  -XI 

N  J  LO  ij 


U+XV  -(U+XV)X  +  (M+XN) 
V  -VX+N 


fU+XV  X(N-VX)  +  M  -  UXl  (rearrange  l-2th 
L  V  N  -  VX  J  entry)  (4.28) 


'  f » (A) 
* 


(from  (4.27)) 


(4.29) 


the  2-2th  entries  of  (4.28)  and  (4.29)  imply 

N  -  VX  =  0  (4.30) 

and  for  the  l-2th  entries, 

0  =  X(N-VX)  +  M  -  UX 

=  OX  +  (M  -  UX)  (from  (4.30)) 

=  M  -  UX 

Thus  a  solution  of  (4.15)  will  satisfy  the  pair  of  equations 
(4.25).  ■ 


The  example  -for  Theorem  4.10  will  serve  to  illustrate 


Theorem  4.11s 


Example  of  Theorem  4.11.  From  (4.22)  begin  by  finding 


the 

characteristic  equation 

for 

r  z 

5  1 

-B  = 

L  z+s 

z-sJ 

f-B 

(p)  =  det (-B  -  pi)  =  det 

rz-p 

s 

lz+s 

z-s 

=  (z-u)  (z-s-p>  -  s(z+s> 

■  z3  -  zs  -  zp  -  J)Z  +  Ps  +  p3  -  sz  -  s2 
=  (z3  -  s3  -  2sz )  +  <s  -  2z)p  +  p3 


Using  the  program  in  Appendix  A,  f-B(R)  = 


-2z3+3sz 

0 

0 


3s3-2sz 

4z=-6sz+s: 

0 

0 


— 2sz3+2s3z+3s3 
-5s3z  +2z  3+s3+sz  3 


0 

0 


5s3z-2s3-2sz3 
- 1 2sz  2+4^3+ 1 Os3z -s3 


[U  Ml 

Lv  Nj 


The  pair  of  equations  which  need  to  be  satisfied  are 


N 

M 


-  VX  ■ 

-  UX  * 


0 

0 


(4.31) 

(4.32) 


(4.31)  doesn't  contribute  to  the  solution,  and  so  the 
solution  to  (4.15)  will  be  found  in  the  solution  space  of 

(4.32)  M  =  UX  or 

-2sz3+2s3z+3s3  5s3z-2s3-2sz3  i 

-5s 2 z +2z  3+s3+sz  3  - 1 2sz 3+4z  3+ 1 0s=z -s3J 

=*  (■  s2  3s3-2sz  I  rw  xi 

[-2z3+3sz  4z3-6sz+s3J  [y  zj  ■ 

To  summarize,  given  the  Lypanov  equation 

AX  -  XB  -  C  A ,  B ,  C ,  X  t  <t"Xr’(D 


and  fA(p),  f„(p)  the  characteristic  equations  of  A  and  B 
respectively,  then  X  will  be  a  common  solution  to  the  pair 
of  equations 

M  =  -XN 
M'  =  U'X 


where 


r  u 

Ml 

and 

f  b(R)  =  fU' 

M'l 

Lv 

NJ 

LV' 

N'J 

Riccati  Equation.  Attention  will  now  be  turned  towards 
(4.2),  where  A,B,C,D,X  f  tr*Xr'(p)  As  with  the  Lypanov 
equation,  a  few  insights  (not  necessarily  connected  with 
solving  the  Riccati  equation)  will  serve  to  motivate 
the  techniques  used  to  solve  (4.2). 

If  X  is  an  n  X  n  matrix  and  I  =  I^x^,  then  the  matrix 

E .[;  ‘] 


always  has  a  determinant  equal  to  -1  (this  follows  from  the 
permutational  definition  of  the  determinant) ,  and  is  thus 
uniutodular .  E's  unimodular  inverse  is 


E' 


since  ee'  -  [i!  a  [?  .a  ■  [i  ?] 

and  E'E  *  I  likewise  follows. 


Given  the  matrix 


R  =  j 

where  A,B,C,D  are  nXn  matrices, 
ERE'1  *  fX  II  TA  B1  fO 


Lc  J 


(4.34) 


X  l] 

r  A  B1  r  0 

I 

I  OJ 

[c  DJ  l  I  - 

■X 

XA+C 

XB+Dl  r  0 

I 

A 

B  J  [  I  - 

■X 

XB+D 

XA+C— XBX— DX  "J 

B 

A-BX  J 

XB+D 

— XBX-DX+XA+C ' 

B 

A-BX 

=  fXB+D  -XBX-DX+XA+C1  (4.35) 

L  B  A-BX  J 

It  can  be  seen  that  the  l-2th  entry  of  ERE_X  is  in  the  -form 
of  the  Riccati  equation  (4.2).  However,  to  get  the  l-2th 
entry  to  agree  exactly  with  (4.2)  which  is 

XDX  +  AX  +  XB  +  C  -  0 

*  -XDX  -  AX  -  XB  -  C  (4.36) 

it  will  be  necessary  to  modify  the  R  matrix  (4.34).  The 
modification  is  implemented  by  comparing  the  form  of  the 
l-2th  entry  of  (4.35) 

— XBX  -  DX  +  XA  +  C  (from  R) 
with  the  the  second  equation  in  (4.36) 

-XDX  -  AX  -  XB  -  C  (original) 

By  matching  the  coefficient  matrices  A,B,C,D  for  like 
expressions  involving  X,  the  suggested  change  is  from 


-r  B] 


(4.37) 


So  with  this  new  form  of  R, 
ERE~l  =  rx  n  r  -b  d  1  ro 


rx  ii 

-B  D  ’ 

ro  1 1 

LI  0. 

.  -C  A  . 

I  -xj 

XD+A 

D 


-XB-C-XDX-AX 
-B  -  DX 


(4.38) 


If  X  is  a  solution  of  the  l-2th  entry,  then 

-XB  -  C  -  XDX  -  AX  =  0  =  XDX  +  AX  +  XB  +  C 
and  X  is  a  solution  of  (4.2). 

The  next  insight  stems  from  the  procedure  followed  in 
solving  the  Lypanov  equation 

AX  -  XB  =  C 

In  Theorems  4.10  and  4.11,  a  solution  to  the  Lypanov 
equation  was  found  among  the  solutions  to  the  pair  of 
equations  (4.17) 

U  +  XV  -  0 
M  +  XN  *  0 

(derived  from  using  the  characteristic  equation  of  A  f«)  and 
(4.25) 

N '  -  V'X  =  0 
M '  -  U'X  =  0 

(derived  from  using  the  characteri stic  equation  of  B  f ») . 
Although  only  one  equation  from  each  pair  turned  out  to  be 
significant  in  solving  (4.1),  the  original  pairs  point  to 
the  solution  of  the  Riccati  equation.  To  obtain  (4. 17) , 
notice  that  if 


ri 

XI 

ru 

Mi 

S 

ru  +  xv 

M  + 

lo 

ii 

lv 

NJ 

l  V 

N 

f°  °1 

l  V  N  J 


.*•  «.*’•  »"•  .  * 
*  V-  X." 


then  matching  the  top  rows  will  yield  <4. 17).  The  matrix 


[l  J] 

was  picked  because  it  helped  generate  the  Lypanov  equation 
in  (4.10).  As  before, 


f«(R) 

S 

ru 

Ml 

,  R  = 

FA 

Cl 

Lv 

nJ 

[o 

bJ 

The  development  of  the  pair  (4.25)  is  similar:  set 

r u  mi  ri  -xi  =  r  u  -ux  +  m] 

IV  N  J  LO  I  J  l  V  -VX  +  NJ 


equal  to 


ru  o' 

|  where 

f»(R>  = 

ru 

M  I 

Lv  o. 

L  v 

N  J 

to  obtain  the  pair 

(4.25) . 

The  same  approach  will  now  be  used  to  generate 
solutions  to  the  Riccati  equation.  Let  f(x)  be  any 
polynomial  in  the  indeterminate  x  with  degree  i  1. 

From  (4.37) 

R  *  f  -B  D  ] 

L  -C  A  J 

Let  f  (R)  »  fU  Ml  (4.39) 

LV  Nj 

Now  the  crucial  step:  as  if  developing  the  pair  (4.17),  let 
X  be  such  that 

Ef  (R)  =  fX  I]  ru  Ml 

1 1  oJ  [v  nJ 

=  rxu+v  XM+N1  =  ro  0]  =  s 

[  U  M  J  [u  MJ  (4.40) 

That  such  an  X  exists  is  motivated  by  the  previous 

discussion  concerning  the  Lypanov  equation.  From  (4.40) 


f  (R)  =  E-1S  or 


f  (R)  » 

fU 

mi  =  ro 

n  r  o 

°l  =  r 

U  M  I 

(4.41) 

lv 

Nj  Li 

-xj  lu 

M 1  L- 

-XU  -XMJ 

Because 

f  (x> 

is  a  polynomial  in  x 

,  xf  (x )  = 

f (x ) x  and 

so  from 

(4.41) 

Rf  (R)  = 

r-B 

D 1  r  u 

M  1  = 

r  -BU-DXU 

-BM-DXMl 

(4.42) 

l-c 

AJ  L-xu 

-xmJ 

[  — CU— AXU 

-cm-axmJ 

and 

f  (R)  R  = 

r  u 

M  1  r-B 

D1  *  r 

-UB-MC 

UD+MA  ] 

(4.43) 

l-xu 

-XM  J  L-c 

aJ  [ 

XUB+XMC  - 

-XUD— XMA J 

Since  (4.42)  and  (4.43)  are  the  same  matrix 


BU  +  DXU  ■  UB  +  MC  (4.44) 
BM  +  DXM  =  -UD  -  MA  (4.45) 
CU  +  AXU  =  -XUB  -  XMC  (4.46) 
CM  +  AXM  *  XUD  +  XMA  (4.47) 


Within  the  above  -four  equations  lie  two  solutions  to 
the  Riccati  equation  (4.2).  From  (4.44)  and  (4.46), 

CU  +  AXU  =  -X (UB  +  MC)  =  -X (BU  +  DXU)  =  -XBU  -  XDXU 
and  so 

CU  +  AXU  +  XBU  +  XDXU  =  0 

=  (C  +  AX  +  XB  +  XDX ) U  (4.48) 

If  det(U)  is  a  nonzero  complex  number,  then  U_1  exists  and 

from  (4.48) 

XDX  +  AX  +  XB  +  C  ■  O 

and  X  satisfies  the  Riccati  equation  (4.2).  The  second 
solution  uses  (4.45)  and  (4.47): 

CM  +  AXM  -  X(UD  +  MA)  =  X(-BM  -  DXM)  =  -XBM  -  XDXM 
and  so 

CM  +  AXM  +  XBM  +  XDXM  =  0  =  (C  +  AX  +  XB  +  XDX)M 


4.20 


If  det (M)  is  a  nonzero  complex  number,  M_1  exists  and  so 


XDX  +  AX  +  XB  +  C  -  0 
another  solution  to  (4.2). 

The  following  theorem  due  to  Jones  has  now  been 
proven* 

Theorem  4.  12*  Given  A,B,C,D  f  4"Xr,(r)  and  f(x)  a 
polynomial  with  degree  >  1,  let 


r _B 

D1 

and 

f  (R)  = 

ru 

Ml 

L-c 

aJ 

lv 

NJ 

If  X  i  4r,x"(r)  satisfies 


XU  +  V  ■  0 
XM  +  N  =  0 


(4.49) 


and  det(li)  it-  <0> ,  then  XDX  +  AX  +  XB  +  C  =  0.  Also, 
if  X  (  4"XrMD  satisfies  the  pair  (4.49)  and 
det  (U)  6  <t  -  f0>  then  XDX  +  AX  +  XB  +  C  =  0. 

Proof s  previous  discussion.  ■ 

A  third  solution  exists  for  (4.2)  but  requires  a 


different 

setup 

akin 

to  developing 

the 

pair  (4.25).  Let 

R  be 

as  in 

i  (4.37)  and 

f (R)  as  in  (4.39) 

.  Let  X  satisfy 

f  (R) 

E  = 

ru 

Ml  r 

0  M  =  f 

« 

U— MX  I 

lv 

NJ  [ 

i  -xJ  L 

N 

V— NX  J 

3 

r M 

°1 

(4.50) 

IN 

oj 

That 

such 

an  X 

exists 

is  motivated 

by  the  development 

of  the 

pair 

(4.25)  for 

the  Lypanov  equation. 

(4.50)  implies 

f  (R) 

3 

ru 

Ml  * 

[!  Sit? 

I] 

=  [MX  Ml 

(4.51) 

lv 

N  J 

oj 

LNX  NJ 

Since  Rf (R)  =  f (R) R 


Rf  (R)  =  r-B  D]  [MX  Ml 

[-C  A  J  L  NX  Nj 

-  f-BMX  +  DNX  -BM  +  DN  1  (4.52) 

L-CMX  +  ANX  -CM  AN  J 

and 


f  (R)  R 


MX  Ml  r-B 
.NX  Nj  [-C 

f-MXB  -  MC 
L-NXB  ~  NC 


S] 

MXD  +  MA  1 
NXD  +  NAJ 


Equating  last  rows  o-f  (4.52)  and  (4.53) 


(4.53) 


NXB  +  NC  ■  CMX  -  ANX  (4.54) 

NXD  +  NA  =  -CM  +  AN  (4.55) 


It  -follows  that 

NXB  +  NC  -  (CM  -  AN) X  =  (-NXD  -  NA)X 
■  -NXDX  -  NAX 


and  so 


NXDX  +■  NAX  -*•  NXB  +  NC  *  0  ■  N(XDX  +  AX  +  XB  +  C) 

If  det(N)  is  a  nonzero  complex  number,  N- 1  exists  and  so 

XDX  +  AX  +  XB  +  C  -  O 


and  a  third  solution  to  (4.2)  has  been  -found  and  thus 

Theorem  4. 13:  Given  A,B,C,D  t  «"x"(r)  and  f(x)  a 
polynomial  with  degree  2  1,  let 


r  -b 

D1 

and 

f  (R)  * 

ru 

Ml 

l  -c 

aJ 

Lv 

Nj 

I-f  X  t  tnXn ( n )  satisfies 


U  -  MX  -  0  (4.56) 

V  -  NX  *  0 

and  det(N)  t  <t  -  <0> ,  then  XDX  ♦  AX  +  XB  +  C  =  0. 


Proofs  previous  discussion.  ■ 

The  next  theorem  presents  a  companion  Riccati  equation 


•  .*  1  ‘  -  '  V  * 


4.22 


to  that  of  (4.2)  given  a  particular  polynomial  -f(x) 


Theorem  4.  14;  Let  f(x)  be  a  polynomial  of  degree  >  1 

and 

XDX  +  AX  +  XB  +  C  =  0  A,B,C,D,X  *  <t"x"(p) 
with  R  and  f  (R)  given  by  (4.37)  and  (4.39).  Then 

XMX  +  NX  -  XU  -  V  ■  0 
Proofi  It  has  already  been  shown  that 

e  -  {x  «o]  «.  e-  -  [o 

are  both  uni  modular  and  inverses  of  one  another.  From 

(4.38)  and  the  given, 

ERE"1  *  fXD  +  A  0  1 

[  D  -B  -  DXj 

From  Theorem  4.8, 

f  (ERE-1 )  *  ff(XD+A)  0  1  (4.57) 

L  *  f(-B-DX)J 

=  Ef(R)E-1  (from  Theorem  4.6) 


=  rxu  +  v  xm  +  n]  ro  it 

L  u  m  J  Li  — x  J 

■  [XM  +  N  XU  +  V  -  XMX  -  NX]  (4.58) 

L  M  U  -  MX  J 

Equating  the  l-2th  entries  of  (4.57)  and  (4.58), 

XU  +  V  -  XMX  -  NX  =  0  =  XMX  +  NX  -  XU  -  V  ■ 


SI 


XAXBX 


This  chapter  will  discuss  some  o-f  the  approaches  which 
may  yield  solutions  to  the  third  order  Riccati  equation 


XAXBX  +  XCX  *  XD  +  EX  +  F 


(5.1) 


where  A,B,C,D,E,F,X  (  4"x"(p>,  r  the  tuple  of  variates.  The 
approach  for  (5.1)  will  be  less  general  than  that  for  the 
second  order  Riccati  and  Lypanov  equations  because  of  the 
form  of  (5.1).  Difficulties  stemming  from  this  form  are 
discussed  in  the  n._«  c  chapter. 

□ne  approach  to  solving  (5.1)  is  akin  to  that  taken  for 
the  second  order  Riccati.  Let  A,B,C,D,E,F  be  known  matrices 


in  G''**"  ( r)  and 


R  * 


[A  B] 
1C  DJ 


f  (R) 


fU  Ml 

L  v  nJ 


where  f(x)  is  a  given  polynomial  with  degree  >1.  As  seen 


previously,  if  X  t  $r’*r,(r)f 

[o  3  [U  "]  [l  -3  ■ 


ru  +  xv  -(u 

L  V 


+  XV)  X  +  (XN  +  M)1 
N  -  VX  J 


Suppose  U  +  XV  *  0  (motivated  by  its  appearance  in 
solving  for  the  Lypanov  and  second  order  Riccati  equations). 
Thus 


(5.2) 


Because  f (R)  is  a  polynomial  in  R,  R  and  f (R)  commute.  From 


f <R>  =  r -xv  M 


j -XV  Ml 
L  V  NJ 


A  B]  r-xv 
C  DJ  L  V 


n]  ’  [ 


BV  -  AXV  AM 
DV  -  CXV  CM 


+  BN] 

+  dnJ 


r-xv 

M] 

fA 

B]  * 

r  mc  - 

XVA 

MD  -  XVB] 

L  v 

nJ 

LC 

dJ 

Lva  + 

NC 

VB  +  ND  J 

Equating  the  -first  columns, 


BV  -  AXV  =  MC  -  XVA 
DV  -  CXV  =  VA  +  NC 


(5.3) 

(5.4) 


VA  appears  in  the  right  hand  sides  o-f  (5.3)  and  (5.4).  If 
MC  in  (5.3)  involved  XNC,  then  a  substitution  from  (5.4) 
could  be  made  in  (5.3).  Accordingly,  let 

MC  -  -XNC  -  Z  (5.5) 

Then  (5.3)  becomes 

BV  -  AXV  »  (-XNC  -  E)  -  XVA 
»  — X (NC  +  VA)  -  Z 

*  -X (DV  -  CXV)  -  Z  (from  (5.4)) 

-  -XDV  +  XCXV  -  E 
or 

XCXV  -  XDV  +  AXV  -  BV  -  Z  -  0  (5.6) 

Because  all  terms  except  Z  end  in  V,  let 

E  =  eiV 

(5.6)  then  becomes 

(XCX  -  XD  +  AX  -  B  -  ex ) V  =  0 
If  det(V)  is  a  complex,  non-zero  number  (easily  verified 
since  f (R)  is  known),  then  the  inverse  of  V  exists  and  so 

XCX  -  XD  +  AX  -  B  -  e,  =  0  (5.7) 

Since  solutions  of  the  third  order  Riccati  equation  are 

5.2 


sought,  let  e. 


-XEXFX 


(5.7)  then  becomes 


XEXFX  +  XCX  +  AX  -  XD  -  B  *  0  (5.8) 

and  X  is  a  solution  to  this  third  order  Riccati  equation. 
However,  since  et  =  -XEXFX  and  2  =  exV,  it  follows  that 
2  *  -XEXFXV  and  from  (5.2),  £  =  XEXFU.  From  (5.5)  it 
follows  that 

MC  ■  — XNC  -  L  -  — XNC  -  XEXFU 

To  summarize,  if 


r  a 

B1 

9 

f  (R)  * 

Fu 

Ml 

l  c 

dJ 

L  v 

N  J 

det(V)  a  non-zero  complex  number  and  X  satisfies  the  pair  of 
equations 

U  ♦  XV  =0 

-XEXFU  -  XNC  «  MC  (5.9) 

9  then  X  is  a  solution  to  the  third  order  Riccati  equation 

XEXFX  +  XCX  +  AX  -  XD  -  B  -  0  (5.10) 


If  det(FU)  is  a  non-zero  complex  number,  then  (5.9)  can  be 
reduced  to  a  second  order  Riccati  equation  of  the  form 
solved  in  the  previous  chapter. 


Another 

sol 

ution  to 

(5 

.  10) 

can  be  gleened 

using 

the 

above  approach. 

Let  XN  •* 

■  M 

a  0 

.  Then 

Rf  (R)  = 

rA 

b]  ru  - 

■XN 

fAU  + 

BV 

BN  - 

AXNl 

(5. 

11) 

lc 

dJ  lv 

N  , 

J 

DV 

DN  - 

CXNJ 

and 

f  (R)  R  = 

ru 

-XNl  [A 

B1 

— 

fUA  - 

•  XNC 

UB 

-  XND” 

|  (5. 

12) 

lV 

N  1  Lc 

DJ 

IvA  + 

■  NC 

VB 

+  ND  J 

Equating  the 

second  columns 

of 

(5.11) 

and 

(5. 12)  , 

-AXN 

+ 

BN  = 

UB  - 

XND 

(5. 

13) 

-CXN 

+ 

DN  = 

VB  + 

ND 

(5. 

14) 

Noticing  that  ND  appears  in  the  right  hand  sides,  let 


UB  =  -XVB  +  L 


(5.  15) 


Then  <5. 13)  becomes 

-AXN  +  BN  ■  (-XVB  +  £)  -  XND 
=  -X(VB  +  ND)  +  £ 

=  — X (DN  -  CXN)  +  E 

*  -XDN  +  XCXN  +  E  (5. 16) 

and  so 

XCXN  -  XDN  +  AXN  -  BN  +  E  =  O  (5. 17) 

Let  E  *  exN.  Then  (5.17)  becomes 

(XCX  -  XD  +  AX  -  B  +  ei)N  =  O 
If  det(N)  is  a  non-zero  complex  number,  then  N  has  an 
inverse,  and  so 

XCX  -  XD  +  AX  -B+ei  =0  (5.18) 

If  ex  =  XEXFX  then 

XCX  -  XD  +  AX  -  B  +  XEXFX  *  0 
and  another  solution  to  the  third  order  Riccati  has  been 


found.  Since  E  =*  e,N,  E  =  XEXFXN  =  -XEXFM.  (5.15)  then 
becomes 


UB  =  -XVB  -  XEXFM 


Thus,  if  X  is  a  common  solution  to  the  pair  of  equations 

0  -  XN  +  M 
UB  =  -XVB  -  XEXFM 


and  det(N)  is  non-zero  complex,  then 

XEXFX  +  XCX  +  AX  -  XD  -  B  -  0 
Another  approach  to  solving  (5.1)  extends  the  solution 
via  similarity  transf ormat i on  used  in  the  matrix  equations 
presented  herein.  Though  this  approach  fails  to  generate 


sets  of  equations  (whose  common  solutions  satisfy  (5.1)),  ?t 
does  present  relations  among  the  coefficient  matrices  for  a 
restricted  case  of  (5.1). 

Let  A,B,C,D,E,F,G,H  i  $r,Xr,(r).  Then  given  X  t  <tr’Xr*(r), 

v  ?]  [i  n  fs  j]  is  a  v  -a 

=  fE+XG  F+XH1  TA+XC  -AX-XCX+B+XD] 

L  G  H  i  [  C  D-CX  j 

=  3 

By  multiplying  these  matrices,  the  following  entries  of  3 
are  obtained: 

1-1:  EA  +  EXC  +  XGA  +  XGXC  +  FC  +  XHC  (5.19) 

1-2:  — XGXCX  -  X (GA+HC) X  -  (EA+FC) X  +  X (GB+HD)  + 

(EB+FD)  +  (EXD— EXCX+XGXD)  (5.20) 
1-3:  GA  +  GXC  +  HC  (5.21) 

1-4:  GB  +  GXD  -  GAX  -  GXCX  +  HD  -  HCX 

Of  the  entries,  the  only  one  in  which  the  "cubic*'  of  the 

third  order  Riccati  appears  is  (5.20),  that  is,  XGXCX.  If 

it  weren't  for  the  last  term  in  paranthesis  in  (5.20),  the 

form  of  (5.20)  would  be  the  same  as  that  of  (5.1).  In  order 

to  eliminate  this  last  term  and  preserve  the  structural 

integrity  of  (5.20),  let  D  *  E  *  0.  It  follows  that 

5  =  fXGA+XGXC+FC+XHC  -XGXCX-X (GA+HC) X-FCX+XGB]  (5.22) 

L  GA+GXC+HC  GB-GAX-GXCX-HCX  J 

If  X  were  a  solution  to  3(1,2),  then  X  would  be  a  solution 

to  a  third  order  Riccati  equation  of  the  form 

XAXBX  +  XCX  +  XD  +  EX  -  0  (5.23) 

which  is  a  special  case  of  (5.1).  From  the  remarks 

preceeding  (5.19),  it  follows  that  (since  D  =  E  =  0) 


■3"  'V 


v« 


5  = 


u  a  [°s  a  [i  a  c  si  [o  -a 

=  f XGA+XGXC+FC+XHC  0  1 

L  GA+GXC+HC  GB— GAX-GXCX-HCX  J 

Multiplying  the  three  middle  matrices  in  (5.24), 


[S  h]  [l  J]  [cfl  S]  ■  [ 


FC 

GA 


HC 


GXC 


0 

GB 


(5.24) 


(5.25) 


j  (5.26) 


I-f  -f(x)  is  any  polynomial  o-f  degree  £  1,  then  -from  (5.24) 
and  Theorem  4.6, 


f  <§>  = 

[1  a  f(5-’  ll 

X  FH 

1 

and  from  Theorem  4.8 

f  <§)  =  ri 

XI  f f (FC)  0  1 

Ij  L  *  f  (GB)  J 

ri  -xi 

(5.27) 

Lo. 

Lo  ij 

r 


XGA+XGXC+FC+XHC) 


f  (GB-GAX-GXCX-HCX ). 


Since  f(J)  and  f(3x)  are  similar  matrices,  they  have  the 
same  determinants.  From  the  permutational  definition  of  the 
determinant  (2.13), 

det (f (J) )  *  detCf (XGA+XGXC+FC+XHC) 3  * 

detCf (GB-GAX-GXCX-HCX) 3  (5.28) 


and 

det (f ( Ji ) )  =  detCf (FC) 1  *  detCf (GB) 3  (5.29) 

and  thus 

detCf (FC)1  *  detCf (GB) 3  =  detCf (XGA+XGXC+FC+XHC) 3  * 

detCf (GB-GAX-GXCX-HCX) 3 

The  following  theorem  has  now  been  proven: 

Theorem  5.1:  Given  A,B,C,F,G,H  (  $ "XrMr>  and 

X  f  t,-,Xr’(|-')  a  solution  to 


5.6 


-XGXCX-X (GA+HC) X— FCX+XGB  *  0 

then 

det  Ef  <FC)  1  *  detEf  (GB)  1  -  det Cf ( XGA+XGXC+FC+XHC) 3  * 

detEf  (GB-GAX-GXCX— HCX )  3 

Proof :  previous  discussion.  ■ 

Theorem  5.1  may  be  extended  to  a  Riccati  equation  of 
the  form 

XAxXAaX  +  XA3X  +  A*X  +  XAo  =  0  (5.30) 

To  get  this  equation  in  the  form  (of  the  given  in 
Theorem  5.1) 

-XGXCX-X (GA+HC) X-FCX+XGB  =  0 
equating  coefficient  matrices  is  done.  Thus 


Ax  =  -G  (5.31) 
A2  =  C  (5.32) 
A3  =  -(GA  +  HC)  =  AiA  -  HA2  (5.33) 
A*  =  -FC  =  -FA3  (5.34) 
A=  ■  GB  -  -AiB  (5.35) 


Since  G  and  C  are  readily  known,  it  remains  to  determine 
A,H,F  and  B.  From  the  above  equations,  one  may  have 
considerable  latitude  in  choosing  A,H,F,B.  From  (5.35),  B 
satisfies  the  equation  A*B  =  -A=.  Likewise,  F  satisfies 
the  equation  FA=  =  -A*.  Lastly,  A  and  H  satisfy  the  Lypanov 
equation  (5.33) 

A3  =  AxA  -  HA3 

By  proper  choice  of  A,B,F,H,  the  challenge  of  finding 
solutions  to  (5.30)  via  Theorem  5.1  may  be  made  easier. 


inclusions  and  Recommendations 


This  paper  sought  the  solution  X  to  -four  types  of 
matrix  equations:  the  linear  equation 

AX  =  b 


the  Lypanov  equation 


AX  -  XB  *  C 

the  second-order  Riccati  equation 

XDX  +  AX  +  XB  +  C  ■  0 
and  the  third-order  Riccati  equation 

XAXBX  +  XCX  +  XD  +  EX  +  F  =  0 
where  the  entries  of  all  matrices  are  restricted  to  being 
multivariate  polynomials  over  the  complex  numbers.  The 
method  of  solution  involved  two  phases: 

1.  identify  a  similarity  transf ormation  on  a  matrix  £ 
in  which  is  embedded  the  coefficient  matrices  of 
the  above  equations.  The  transf ormation  gives  a 
matrix  whose  entries  include  the  equation  being 
sol ved . 

2.  identify  a  polynomial  which,  when  its  argument  is 
the  matrix  Z  in  <1>,  gives  a  matrix  whose  entries 
yield  pairs  of  linear  equations.  The  common 
solutions  of  these  pairs  of  equations  will  contain 
those  of  the  given  matrix  equation.  The  choice  of 
polynomial  may  yield  different  solutions. 

Since  the  search  for  solutions  ends  in  solving  a  linear 
equation,  a  method  was  presented  which  may  identify  the 
general  solution  of  the  given  linear  equation.  The  method's 
success  depends  on  if  the  matrix  A  in  AX  =  b  is  reducible  to 
a  Smith  Form.  The  similarity  transf ormati on  of  (1)  is  given 
by  E  E  E_1  where  E  equals 


(which  is  uni modular)  and  £  for  the  Lypanov  equation  equals 


and  £  for  the  second  order  Riccati  equals 

[:?  2] 

The  third  order  Riccati  presented  difficulties  because  its 
form  wasn't  adaptable  to  the  method  of  similarity 
transformations.  Nevertheless,  a  transformation  could  be 
done  which  might  yield  some  solutions.  When  successful,  the 
solution  to  the  third  order  Riccati  was  a  common  solution  to 
a  linear  and  a  second  order  Riccati  equation.  The 
polynomial  in  (2)  can  assume  any  form  and  when  its  argument 
is  £, 

f  (£)  =  ru  Ml  (6.1) 

L  V  h*  J 

The  submatrices  li,M,V,N  yield  pairs  of  equations  which 
depend  on  the  matrix  equation  being  solved.  These  pairs  in 
turn  yield  solutions  of  the  original  matrix  equation. 

Recommendations.  The  biggest  challenge  is  also  the 
central  concern:  develop  the  mathematical  theory  which  will 
identify  whether  or  not  a  given  matrix,  whose  entries  are 
multivariate  polynomials,  is  reducible  to  Smith  Form. 

Though  the  answer  for  the  two-variate  polynomial  case  is 


known,  the  three  or  more  variate  case  remains  an  open 
question.  It's  felt  that  this  challenge  will  involve 


creative  work  involving  number  theory,  abstract  algebra,  and 
matrix  theory.  Qnce  the  theory  identifying  the  conditions 
under  which  reducibility  to  Smith  Form  is  known,  the  next 
hurdle  will  be  to  develop  a  computer  algorithm  performing 
the  reductions  <the  Smith  Form  may  not  be  unique).  The 
papers  C83  and  Z201  give  approaches  to  the  two-variate  case, 
and  the  papers  C23D  and  C333  give  insights  dealing  with  the 
"computer  programming  algebra”  of  matrix  forms  and 
polynomials.  Perhaps  there  is  an  underlying  linear 
programming  or  network  formulation  to  the  Smith  Form 
reduction. 

In  solving  the  Lypanov  and  Riccati  equations,  it  was 
necessary  to  assume  the  existence  of  the  inverse  to  one  or 
more  of  the  submatrices  in  (6.1).  Of  course,  (6.1)  depends 
on  the  polynomial  f(x)  chosen,  which  in  turn  dictates  the 
types  of  solutions  found  to  the  original  matrix  equation. 
Research  needs  to  be  done  into  the  types  of  polynomials 
yielding  one  or  more  invertible  submatrices  of  (6.1). 

Perhaps  there  are  equivalence  classes  of  polynomials.  Or 
families  of  polynomials.  liaybe  there  is  no  structure.  Or 
perhaps  there  is  a  "minimum"  polynomial  which  (in  some  way) 
generates  all  polynomials  yielding  invertible  submatrices. 
Again,  a  knowledge  of  abstract  algebra  may  prove  valuable: 
much  literature  dealing  with  polynomial  structures  is 
couched  in  abstract  algebraic  terms. 

Much  effort  was  spent  in  finding  the  "proper  embedding" 


■for  the  third  order  Riccati  equation.  Disappointingly 
limited  success  was  realized.  To  see  what  is  meant  by 
"proper  embedding",  notice  that  in 

r  I  XT  [A  B]  [I  -Xl  =  fA  +  XC  -AX  -  XCX  +  B  +  XD] 

[o  IJ  [c  DJ  lo  IJ  [  C  -CX  +  D  J 

a  second  order  Riccati  equation  is  found  in  the  l-2th 

position  of  the  matrix  on  the  right.  If  C  =  0,  then  this 

equation  becomes  the  Lypanov.  Thus  the  Lypanov  and  second 

order  Riccati  equations  are  "embedded"  in  the  matrix  on  the 

right,  which  in  turn  comes  from  the  similarity 

transf ormation  on  the  left.  Since  this  approach  has  proven 

so  successful  in  generating  linear  equations  whose  solutions 

are  those  to  the  Lypanov  and  second  order  Riccati  equations, 

it  is  natural  to  extend  the  similarity  transf ormation  to 

handle  the  third  order  Riccati  equation. 

Unfortunately,  this  is  more  easily  said  than  done. 

Since  the  third  order  Riccati  equation 

XAXBX  +  XCX  +  XD  +  EX  +  F  =  0  (6.2) 

has  six  coefficient  matrices,  the  above  approach  would  have 
to  be  modified  to  handle  these  six.  As  it  stands,  it  is 
geared  to  the  four  coefficients  of  the  second  order  Riccati 
equation.  Is  it  possible  to  work  with  matrices  which  are 
two  sub-matrices  deep  and  wide,  or  is  it  necessary  to  go  to 
higher  order  matrices,  perhaps  three  submatrices  in 
dimension.  If  so,  what  would  the  unimodular  matrix  be  that 
effects  the  similarity  transf ormati on?  Is  the  restriction 
of  unimodularity  unduly  restrictive,  mandatory  though  it 


seems?  Also,  the  first  term  XAXBX  of  (6.2)  presents  a 
challenge:  how  to  introduce  the  middle  X  from  the 

transf ormation  E  £  E-1  while  hopefully  embedding  a  third 
order  Riccati  equation  in  the  resulting  matrix  of 
transition.  Perhaps  a  different  type  of  transf ormation 
needs  to  be  done.  Instead  of  a  single  transform,  a  string 
of  transforms  may  be  the  answer,  for  example 

Et  C1E2S2E3 

where  the  £j  are  matrices  whose  entries  include  the 
coefficient  matrices,  and  the  Ej  contain  the  solution  X.  If 
this  approach  is  indeed  the  way  to  go,  how  are  the  Ej's 
related.  Again,  the  third  order  Riccati  should  be  embedded 
somewhere  in  the  resulting  matrix. 

Or  should  it?  If  not  embedded,  could  other  matrix 
equations  which  are  embedded  as  a  matter  of  course  yield,  in 
some  combination,  some  (or  all!)  solutions  of  the  third 
order  Riccati?  Another  concern  arises:  in  the  approaches 
presented  herein,  a  specially  selected  polynomial  f(x) 
played  the  crucial  role  of  generating  pairs  of  equations 
whose  common  solutions  included  those  of  the  matrix  equation 
in  question.  For  the  third  order  Riccati,  how  is  the 
polynomial  to  be  selected?  How  shrjld  the  pairs  be 
generated,  if  indeed  this  pattern  holds?  Perhaps  triplets 
of  equations  arise  instead  of  pairs. 

The  unsettling  thought  occurs:  maybe  similarity 
transforms  are  not  powerful  enough  for  the  third  order 


Riccati,  and  whole  new  approaches  are  in  order.  After  all, 
current  research  into  these  matrix  equations  is  still  in  the 
pioneering  stages.  Hopefully,  the  successful  resolution  to 
the  third  order  Riccati  equation  will  point  the  way  to 
solving  higher  order  Riccati  equations.  However,  a  spectre 
appears:  in  solving  for  the  roots  of  a  polynomial  in  one 

variable,  the  insolvability  of  the  quintic  was  demonstrated 
by  Abel  in  the  early  nineteenth  century.  Could  there  be  a 
similar  obstacle  ahead  for  matrix  equations?  In  any  case, 
finding  solutions  to  matrix  equations  will  prove  challenging 
and  endlessly  fascinating — as  well  as  having  immediate 
practical  applications. 

The  last  recommendation  is  a  minor  one:  create  a  new 
segment  (in  the  computer  program  of  Appendix  A)  which  will 
compute  the  determinant  of  a  given  square  matrix  having 
multivariate  polynomial  entries.  The  method  would  have  to 
keep  calculations  down  to  a  minimum,  since  finding 
determinants  is  computationally  intensive.  One  approach, 
which  seems  to  be  the  quickest  way  to  find  a  matrix 's 
determinant,  was  published  by  G.  Macloskie  C223  in  1904. 

This  method  could  be  adapted  to  that  of  finding  determinants 
for  matrices  with  polynomial  entries. 


6.6 


Appendix  A: 


~ming  matrix 


f 


The  following  BASIC  program  performs  operations  on 

rectangular  matrices  having  multivariate  polynomial  entries 

over  the  real  numbers.  Via  menu  options,  the  interactive 

program  allows  the  user  to 

-create  a  matrix 
-view  a  matrix 
-transpose  a  matrix 
-add  two  matrices 
-multiply  two  matrices 

-extract  the  U,M,V,N  matrices  from  a  given  matrix 
The  size  of  matrices  which  the  program  can  handle  is  limited 
only  by  the  amount  of  memory  available  on  the  computer  and 
the  parameters  in  the  DIMENSION  statement.  The  program  was 
designed  on  an  IBM  PC  (DOS  3.0>  and  was  intended  to  assist 
the  thesis  effort. 

The  program's  logic  hinges  upon  the  way  it  recognizes  a 
polynomial  in  n-variates:  as  a  set  of  (n+1) -tuples  each  of 
whose  entries  come  from  the  coefficients  and  exponents  of 
the  polynomial.  To  obtain  the  tuples,  the  entire  polynomial 
is  re-expressed  in  "standard  form"  as  a  sum  of  terms,  each 
having  all  variate  and  exponent  positions  appear.  For 


instance,  the  polynomial  in  the  variates  x  and  y 

2xay  —  17xya  -  3y  +  7 


is  re-expressed  in  standard  form  as 


(A.  1) 


2x=y1  +  <“17x  lya)  +  (-3x°yM  +  (7x°y«)  (A. 2) 

and  is  the  sum  of  the  four  terms  2x=yl,  -17x1ya,  -3x°y1, 

7x°y°.  Though  the  order  of  the  variates  <i.e.,  x  before  y 


.  a V-V-V- 


'll, 


versus  y  before  x)  is  arbitrary,  the  order  must  be 
consistent. 

□nee  in  standard  form,  each  term  is  put  into  its 
tuple  representation 

(coefficient,  exponent  *  ,...  , exponent,-.)  (A. 3) 

with  the  polynomial  represented  by  the  collection  of  the 
tuples.  For  the  above  example: 


TERM 

3-tuple 

2x2y 1 

<  2,2,1) 

17x 1y= 

(-17,1,5 

■3x°y 1 

(-3,0,1) 

7x°y° 

(  7,0,0) 

and  so  the  polynomial 

2x2y  -  17xya  -  3y  +  7 
is  represented  by  the  set 

<  (2,2,1),  (-17, 1,5),  (-3,0,1),  (7,0,0)  > 
Similarly,  the  polynomial 

42xyuvz  -  2z3yx 

is  re— expressed  in  the  standard  form  (variate  order  is 
x  »y,z ,u,v) 

42x 1y1z1u1vl  +  (— 2x 1y1z2u°v0) 
and  thus  the  set  of  6-tuples 

<  (42,1,1,1,1,1),  (-2, 1,1, 2, 0,0)  > 

Addition  and  multiplication  between  two  polynomials  are 
easily  expressed  in  terms  of  their  tuple  representations 


provided 


i 


-all  tuples  have  the  same  number  o-f  entries 
— coordinate  positions  all  correspond  to  the  same 
variable 

Addition.  Given  two  terms  in  n-variates,  let  (a, a)  and 
(b,fl)  be  their  respective  (n+1) -tuple  representations  where 
a,b  are  the  coefficients  o-f  each  term  and  a,B  the 
exponent  vectors.  For  example,  i-f  the  terms  are  3x2y  and 


7x^y=*  f  then 


3x=y  %  <3,2,1)  ,  a  -  3,  a  =  (2,1) 
7x’y®  %  (7,9,5)  ,  b  >  7,  B  ■  (9,5) 


Tuple  addition  is  defined  by 


(a, a)  #  (b,fi) 


<  <a, a) ,  <b,fi>  > 
(a+b , a) 


a  *  fi 


a  ■  o 


(A. 4) 


The  operation  #  is  clearly  commutative.  To  add  two 

polynomials,  the  tuples  representing  the  terms  are  added 

according  to  <A.4).  For  example: 

<2x2y  -  9xy*)  +  <4xy*  +  xzy) 

*  C  (2,2, 1)  ♦  (-9,1,2)]  «  [(4,1,2)  =if=  (1,2,1)] 


then,  rearranging  terms 

=  C (2,2, 1 )  (1,2,1) ] 

=  (3,2,1)  =11=  (-5,1,2) 

*  f  (3,2,1),  (-5,1,2)  > 

*  3x2y  -  5xy2 


C  (-9 ,1,2)  *  (4,1,2)] 


Multiplication  between  two  tuples  is 


defined  by 


<a,a)*(b,B)  *  <ab,a+B) 


(A. 5) 


where  ff  +  B  is  the  familiar  operation  of  addition  between 
two  vectors  in  Euclidean  n-space.  For  example, 


Using  (A.  5)  the  tuple  representation  of  the  product 
between  two  polynomials  is  quickly  -found.  Consider  the 
product  o-f  the  n-term  polynomial 

TERM  x  x  +  TERMir.  (A.  6) 

and  the  m-term  polynomial 

TERMai  +  ...  +  TERM3m  (A. 7) 

which  is 

TERM* i (TERM3i  +  . . .  +  TERMa™.)  +  ...  (A. 8) 

+  TERMi„ (TERM=i  +  ...  +  TERMa™,) 

*  ( TERM i i TERMa i  +  ...  +  TERM i i TERMam )  +  ... 

+  (TERMi„  TERMai  +  ...  +  TERM  i„  TERMa™,) 

I-f  the  value  o-f  TERMj*  is  now  changed  to  the  tuple  form 

of  the  polynomial  element  TERM;,*,  then  the  product  of  (A.  6) 

and  (A. 7)  expressed  in  tuple  form  is 

=  (TERMi  i  •  TERMai  =11=  .  .  .  #  TERM  1 1  *  TERMa™, )  ^  ... 

=ir  (TERM*,, ‘TERMai  ^  ...  #  TERM  x  „  TERMa™,  > 

EXAMPLE; 

(x2y3  +  12x,*y=>  (x6y7  +  llx'V*) 

~  (x2y3)  (x^y7)  +  (x2y3)  ( 1  lxQy*») 

+  (12x*y=)  (xV)  +  (12x*y=»)  (llx'V7) 

*  C (1 ,2,3) • (1,6,7) 3  #  E (1,2,3) • (11,8,9)1 

nr  C  (12,4,5)  •  (1,6,7)  1  ^  C  (12,4,5)  •  (11,8,9)  1 

=  (1,8,10)  nr  (11,10,12) 

#  (12,10,12)  nr  (132,12,14) 

=  (1,8,10)  nr  (23,10,12)  (132,12,14) 

*  C  (1,8,10),  (23,10,12),  (132,12,14)  > 

;;  x«yio  +  23x10yi=  +  132x13y**  ■ 

Network.  Put  in  a  nutshell,  the  program  keeps  track  of 
tuple  operations  using  a  system  of  pointers  (program  lines 


1150-1490) 


That  is,  a  network  tree  is  constructed  which 


represents  the  results  of  tuple  operations.  This  section 
assumes  knowledge  o-f  networks  as  given  in  C10:91-124D. 


The  fundamental  insight  is  that  a  polynomial  can  be 
represented  by  a  tree.  Specif ically,  the  results  of  tuple 
operations  is  maintained  in  a  network  tree  whose  nodes 
represent  a  particular  variate  in  a  polynomial  term.  Nodal 
"potentials"  are  the  exponent  of  the  corresponding  variate. 
Terminal  nodes  of  the  polynomial 's  tree  have  a  second 
potential  whose  value  is  that  of  the  coefficient  of  a 
particular  term.  The  i ‘th  depth  of  the  tree  corresponds  to 
the  i 'th  variate.  For  example,  the  polynomial 

3x*yz  13  -  7x  l3y"»z  +  8x  1 3y ■’’z 2:2 
is  represented  by  the  tree  in  Figure  1.  All  the  nodes  at 
depth  one  are  x  variates;  nodes  at  depth  two  are  y  variates; 
nodes  at  depth  three  are  z  variates.  The  numbers  in  CD  are 
node  potentials,  and  numbers  in  CD*  are  the  coefficient  of 
the  given  polynomial  term.  From  this  example,  it  can  be 
seen  that  a  polynomial  term  corresponds  to  a  unique  path 
from  node  0  to  a  given  terminal  node  at  depth  three,  e.g., 

the  term  8xl2y"*z==  has  the  nodal  path  0-4-5-7. 

In  performing  a  tuple  operation,  the  order  in  which  the 
branches  are  built  are  term  by  term.  The  following  example 
illustrates  the  sequence  in  which  the  program  builds 
branches. 

EXAMPLE;  Construct  the  tree  representing 

X2(l  -t-  y)  +  y(x2  +  y) 

A. 5 


(A. 9) 


x=(l  +  y)  +  y(x2  +  y)  =  x=  +  x=y  +  yx2  +  y2 

=  X2y°  X=y 1  4"  y 1  x 2  +  x°y= 

Since  there  are  2  variates  x,y  the  tree  will  have  2  depths. 
The  first  term  is  x2y°  which  contributes  the  first  branch 
0-1-2  in  Figure  2.  The  second  term  x2y‘  contributes  the 
branch  0-1-3  in  Figure  3.  Note  that  the  second  potential  of 
node  3  <CJ*)  is  1:  the  coefficient  of  the  second  term  x2yx. 
Because  the  third  term  of  (A. 9)  is  also  x2yx,  it  doesn't 
contribute  a  branch  to  the  tree.  However,  the  second 
potential  at  node  3  is  increased  by  1:  the  coefficient  in 
the  third  term  x2y»  (Figure  4> .  The  fourth  term  x°y2 
contributes  the  branch  0-4-5  in  Figure  5.  Since  there  are 
no  more  terms  in  (A. 9),  the  tree  in  Figure  5  represents  the 
polynomial  x2(l  -►  y)  +  y(x2  +  y) 

Manual  pquu  of  data.  When  creating  a  matrix  via  the 
menu  prompt  tc  ’huild  a  matrix",  polynomials  are  entered 
term  by  term  in  their  tuple  form  with  no  commas  between 
tuple  entries.  Each  tuple  MUST  remain  on  the  same  line,  and 
each  tuple  except  the  last  ends  with  a  carriage  return;  the 
last  tuple  ends  with  a  semicolon  followed  by  a  carriage 
return.  The  program  will  then  prompt  for  the  next 
polynomial  entry  and  the  above  protocol  is  repeated.  A 
sample  session  input ing  the  matrix 

Q  =  [x2  +  3xyz  +  z2  4  +  2x  -  3y]  (A.  10) 

L  z*  +  17y*”s  x2y3z“*  J 

is  given  in  Figure  6. 

Output .  The  program  outputs  a  polynomial  in  its  tuple 


'  O  .  - 

\  w-T  s.  l  U-_". 
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UNCLASSIFIED 


NL 


2  2 

Branch  for  x  -1  2x  y 


Fig.  5.  Branch  for 


********■»■»*•»*■*•»**•»***■»*****■»****•»*■»■*********■»*** 

MENU: 

B  Build  a  new  matrix” 

V  View  a  matrix" 

Transpose  a  matrix” 

+  Add  2  matrices” 

*  Multiply  2  matrices” 

U  Extract  U,M,V,N  matrices 

*#*««**#*#*4Ht***####**##**#*********#**##***-»*-»* 

?  B 

Matrix  name?  B:Q 

ROW, COL  dimension?  2,2 

Number  o-f  polynomial  variates?  3 

Build  a  square  diagonal  matrix  whose  diagonal  entries 
are  equal  Y/N?  N 
***  COEF  EXP1  EXP2  . . .  EXPn  *** 
a  semicolon  ;  represents  end  o-f  polynomial 

(1,1) : 

?  1  2  0  0 
?  3  1  1  1 
?  1  0  0  2; 

(1,2) : 

?  4  0  0  0 
?  2  1  0  0 
?  -3  0  1  0; 

(2,1) : 

?  1  O  O  6 
?  17  0  98  0; 

(2,2) : 

?  1  2  3  4; 


✓  f  // 


************************************************ 

MENU: 

B 

Bui Ida 

new 

matrix " 

V 

View  a  matrix" 

# 

Transpose  a 

matrix " 

+ 

Add  2  matrices" 

* 

Multiply 

2  matrices” 

u 

Extract 

U,M 

,VfN  matrices 

1  ***■»***•»*■»**■»*■»****■**■»•»***■»*****«•***********■»■»**  9 

1  ?  V 

I  Screen  output  Y/N? 

Y 

I  Matrix  to 

view?  B: 

•f  A  (R) 

MATRIX  B: f A (R)  IS 

4  X 

4  and  has 

2  polynomial 

variates 

B: f A (R)  (  1 

,  1): 

B:  f  A  <R)  ( 

3, 

1) : 

0  0  0 

0  0  0 

B:  -f  A  <R)  (  1 

,  2): 

B:  f  A  (R)  ( 

3, 

2) : 

0  0  0 

0  0  0 

B: -f  A (R)  (  1 

,  3): 

B:  -f  A  (R)  ( 

3, 

3) : 

-12  1 

2  2  0 

3  3  0 

-2  1  1 

B:  f  A  <R)  (  1 

,  4): 

B:  -f  A  (R)  ( 

3, 

4) : 

-5  3  0 

-3  2  0 

7  2  1 

2  1  1 

-2  1  2 

B:f A(R) (  2 

.  Ds 

B:  -f  A  (R)  ( 

4, 

1) : 

0  0  0 

0  0  0 

B:  f  A  (R)  (  2 

,  2): 

B:  f A  <R>  ( 

4, 

2) : 

0  0  0 

0  0  0 

B-  -f  A  (R)  (  2 

,  3): 

B:  f A  <R>  ( 

4, 

3)  : 

5  1  2 

-111 

-5  3  0 

-3  2  0 

-2  0  3 

2  0  2 

2  2  1 

Bs-fA(R)  (  2 

,  4): 

B:  -f  A  (R)  ( 

4, 

4)  : 

2  12 

5  2  0 

|  8  3  0 

-4  1  1 

-8  2  1 

NOTE:  the  program  has  been  modified  to  fit  within  the 

format  requirements  of  the  AFIT  Style  Guide,  and  so 
doesn't  represent  the  correct  syntax  of  the  program. 
However,  the  adjustments  are  few  and  are  apparent. 


10 

'CPT  Bruce 

W.  Colletti,  AFIT/EN;  1  August  1985 

20 

'The  input 

file  for  this  program  (which  multiplies 

30 

'matrices  having  multivariate  polynomial  entries) 

'has  the  following 

40 

'structure: 

50 

line  1: 

row,  column  dimension  of  matrix  A; 

# 

number  of  variates 

60 

line  2: 

number  of  terms  in  the  A (1,1) 

0 

polynomial 

70 

line  3: 

tuple  corresponding  to  the  1st  term  of 

A  ( 1 , 1 ) 

80 

'  line  4: 

tuple  corresponding  to  the  2nd  term  of 

A  ( 1 , 1 ) 

90 

line  ns 

tuple  corresponding  to  the  last  term 

of  A  ( 1 , 1 ) 

100 

'  line  n+1: 

number  of  terms  in  the  A (1,2) 

• 

polynomial 

110 

'  line  n+2: 

tuple  corresponding  to  the  1st  term  of 

0 

A ( 1 ,2) 

120 

'  line  n+3: 

tuple  corresponding  to  the  2nd  term  of 

• 

A ( 1 ,2) 

130 

line  m: 

tuple  corresponding  to  the  last  term 

• 

of  A ( 1 ,2) 

140 

'  line  m+1: 

number  of  terms  in  the  A (1,3) 

0 

polynomial 

150 

'  etc: 

matrix  entries  are  entered  row  at  a 

time 

160 

CLS  :  CLEAR 

170 

ROUS  ■  5  : 

COLS  ■  5  :  TERM  -  25  :  VAR  *  5 

180 

DIM  A (ROWS, 

COLS , TERM , VAR ) ,  B ( ROWS , COLS , TERM , VAR ) , 

B*(20) 

190 

DIM  F (90) , POT ( 90 ) ,R(90) ,BACK(90) ,M(90) ,E(VAR) 

200 

DIM  C0EFL ( ROWS , COLS ) , COEFR ( ROWS , COLS ) 

210 

'A  contains 

the  elements  of  the  left  matrix  and  B 

'the  elements  of 

220 

'the  right 

(e. g.  ,  A  +  B,  A*B) : 

230 

'A(I,J,K,L) 

the  exponent  of  the  L'th  variate  of  the 

'K'th  term 

of  the 

240 

0 

I-J'th  polynomial  entry  of  A 

250 

' B ( I , J ,K,L) 

the  exponent  of  the  L'th  variate  of  the 

'K'th  term 

of  the 

260 

0 

I-J'th  polynomial  entry  of  B 

270 

'  BACK ( i ) 

the  node-back  pointer  for  node  i 

280 

' COEFL ( i , j ) 

the  number  of  terms  in  the  polynomial 

located  at 

A(i  ,  j) 

COEFR(i,j)  the  number  of  terms  in  the  polynomial 
located  at 


290 

300 


ro  4-  in  «o  rs  cd  o*  o^Nto<fin'Oi^a)0'O^Nt')«tin'ONa)0‘O^Nro4-inoiscD  o*  o  -« om 


E  (i ) 


B(i , j) 

contains  the  exponent  of  the  i *th 
variable  of  the 

product  between  2  polynomial  terms 
F(i)  forward  pointer  for  node  i 
M(i)  the  "summed  coefficient".  This  is  a 
node  in  the  final 

tree  with  no  forward  pointer.  It  is 
the  coefficient  of 

that  term  of  the  polynomial  having  the 
path  from  given 
terminal  node  to  source  node 
POT(i)  potential  of  node  i;  that  is,  the 
exponent  of  the 

variable  represented  by  the  node 
R(i)  the  right  pointer  for  node  i 
PR I NT  " ************************************* " 

PRINT  "MENU:" 

PRINT  "  B  Build  a  new  matrix" 

PRINT  "  V  View  a  matrix" 

PRINT  "  '  Transpose  a  matrix" 

PRINT  "  +  Add  2  matrices" 

PRINT  "  *  Multiply  2  matrices" 

PRINT  "  U  Extract  U,M,V,N  matrices 

PRINT  "* 


11  fe********^********#*********##******* 11 


INPUT  OPT* 

IF  OPT*  -  "B* 
IF  OPT*  «  "V* 
IF  OPT*  »  "U’ 
IF  OPT*  =  -*■ 
IF  ( (OPT*  <> 


THEN  2410 
THEN  2790 
THEN  3470 
THEN  3130 


IF  ( (OPT*  <>  "*")  AND  (OPT*  <>  "+"))  THEN  420 
INPUT  "Left  matrix";  FILEL* 

INPUT  "Right  matrix";  FILER* 

INPUT  "Output  matrix";  FILEO* 

IF  OPT*  »  «+••  THEN  620 

INPUT  "0t  i  R  in  «LEFT*RIGHT";  ALPHA  :  GOTO  630 

INPUT  "0!,B  (  R  in  a LEFT  +  BRIGHT";  ALPHA, BETA 

FILE*  »  FILEL*  :  GOSUB  950 

R0W1  =  ROW  :  C0L1  =  COL  :  VI  =  V 

FILE*  -  FILER*  :  GOSUB  950 

R0W2  =  ROW  :  C0L2  =  COL  :  V2  a  y 

IF  VI  <>  V2  THEN  1730 

IF  ((OPT*  =  «+«)  AND 

( (R0W1  <>  ROW)  OR  (C0L1  <>  COL)))  THEN  1560 
IF  ((OPT*  »  "*••)  AND  (C0L1  <>  ROW))  THEN  1560 
OPEN  FILEO*  FOR  OUTPUT  AS  #2 
IF  OPT*  =  THEN  GOSUB  2140 
PRINT  #2,  R0W1;  C0L2;  V 
FOR  II  =  1  TO  R0W1 
FOR  12  =  i  TO  C0L2 
PRINT  "(";  II;  " 12;  ">" 

GOSUB  1800 

FOR  13  =  1  TO  C0L1 


T2  =  COEFR (13, 12) 


780 

790 

800 

810 

820 

830 

840 

850 

860 

870 

880 

890 

900 

910 

920 

930 

940 

950 

960 

970 

980 

990 

1000 

1010 

1020 

1030 

1040 

1050 

1060 

1070 

1080 

1090 

1100 

1110 

1120 

1130 

1140 

1150 

1160 

1170 

1180 

1190 

1200 

1210 

1220 

1230 

1240 

1250 

1260 

1270 

1280 

1290 


T1  =  COEFL (11,13)  s 
FOR  14  *  1  TO  T1 
FOR  15  ■  1  TO  T2 

PROD  =  ALPHA  *  A<I1 , 13, 14,0)  * 

B( 13, 12, 15,0) 

FOR  16  =  1  TO  V 

E ( 16)  =  A(I1, 13, 14, 16)  +  B  < 13, 12, 15, 16) 
NEXT  16 
GOSUB  1150 
NEXT  15 
NEXT  14 
NEXT  13 
GOSUB  1870 
NEXT  12 
NEXT  II 
CLOSE  #2 
STOP 

'Subroutine  reads  a  matrix  for  +/*  routines 
OPEN  FILE*  FOR  INPUT  AS  #1 
INPUT  #1,  ROW, COL, V 

IF  ((ROW  >  ROWS)  OR  (COL  >  COLS))  THEN  1620 
FOR  I  -  1  TO  ROW 
FOR  J  ■  1  TO  COL 

INPUT  #1,  TERMS 
IF  TERMS  >  TERM  THEN  1680 
IF  FILE*  =  FILEL*  THEN  COEFL(I,J)  -  TERMS 
IF  FILE*  =  FILER*  THEN  COEFR ( I, J)  *  TERMS 
FOR  K  *  1  TO  TERMS 
FOR  L  *  0  TO  V 

IF  FILE*  *  FILER*  THEN  1080 
INPUT  #1,  A ( I , J ,K,L)  s  GOTO  1090 
INPUT  #1,  B(I,J,K,L) 

NEXT  L 
NEXT  K 
NEXT  J 
NEXT  I 
CLOSE  #1 
RETURN 

'Pointer  ordering  subroutine 
1=1:  CNODE  =  F(0)  :  BK  =  0 
IF  CNODE  =  0  THEN  1370 
BASE  3  BK 

IF  POT (CNODE)  =  E(I)  THEN  1430 
BK  »  CNODE 
CNODE  =  R< CNODE) 

IF  CNODE  <>  0  THEN  1190 
NODES  =  NODES  +  1 
POT  (NODES)  =  Ed) 

R (BK)  =  NODES 
BACK (NODES)  =  BASE 
BK  =  NODES 
FOR  J  =  1+1  TO  V 

NODES  =  NODES  +  1 


1300 

1310 

1320 

1330 

1340 

1350 

1360 

1370 

1380 

1390 

1400 

1410 

1420 

1430 

1440 

1450 

1460 

1470 

1480 
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POT (NODES)  -  E(J) 

BACK (NODES)  *  BK 
F (BK)  =  NODES 
BK  =  NODES 
NEXT  J 

M (BK)  =  M (BK)  +  PROD 
GOTO  1490 
NODES  -  NODES  +  1 
F (NODES)  =  0 
F (BK)  -  NODES 
BACK (NODES)  *  BK 
CNODE  -  NODES 
POT (CNODE)  *  E(I) 

BASE  -  CNODE 

1*1+1  'I  is  the  same  as  the  Depth  of  a  node 
IF  (I  >  V)  THEN  M (CNODE)  =  M (CNODE)  +  PROD  s 
GOTO  1490 
BK  *  CNODE 
CNODE  =  F (CNODE) 

GOTO  1170 
RETURN 

'Error  routine  for  too  many  polynomial  variates 
PR I NT  " **************************************** " 
PRINT  CHR*(7);  "DIM  statement  allows  only  VAR; 
"variables. " ; 

PRINT  "  Modify  DIM  statement"; 

STOP 

'Error  routine  for  ill-dimensioned  matrices 
PR I NT  " **************************************** " 
PRINT  CHR*(7);  "Dimensions  incompatable" 

PRINT  "  ";  FILEL*  ;  "  is  ";  STR*(R0W1);  "  X  "; 

STR* (C0L1) 

PRINT  "  ";  FILER*  ;  "  is  ";  STR*(R0W);  "  X  "; 

STR* ( COL ) 

STOP 

'Error  routine  for  matrices  whose  dimensions  are 
'too  large 

PR I NT  " ***************************************** " 
PRINT  CHR*(7);  "Modify  the  dimension  statement  in 
the  program" 

PRINT  ;  "The  dimensions  of  ";  FILE*;  "is  " ; 

STR* (ROW) ; 

PRINT  "  X  STR* (COL) 

STOP 

'Error  routine  for  a  polynomial  having  too  many 
terms 

PR I NT  " ****************************************** " 
PRINT  CHR*(7);  "There  are  ";  TERMS;  "terms  in 
the  " ; 

PRINT  "polynomial,  whereas  the  DIM  statement 
allocates  " 

PRINT  TERM;  "  terms  in  the  polynomial.  Modify 
DIM  stmt. " 


r  '  J*'  ^  m y  '  v,” r-  :  y  *  -  **  ► 


i»y-: 


PR I NT  " ****************************************** " 
PRINT  CHR*<7);  "Polynomial  variates  in  matrices 
aren't  the  same" 

STOP 

'Error  for  finding  U,M,V,N  matrices 
PRINT  CHR$<7>;  "Source  matrix  ";  FILE1$;  "  doesn't 
have  even"; 

PRINT  "  row  and  column  dimensions" 

STOP 

'Subroutine  to  clear  tree  building  variables 
FOR  I  ■  0  TO  NODES 

F  <  I )  =  0  s  POT  (I)  =0  ;  R  ( I )  =0 
BACK ( I )  *  0  •  M<I>  -  0 
NEXT  I 
NODES  »  0 
RETURN 

'Subroutine  prints  polynomial  term 
S  =  0 

FOR  I  ■  1  TO  NODES  'find  #  of  terms  in  polynomial 
IF  M(I)  <>  0  THEN  S  =  S  +  1 
NEXT  I 

IF  S  ■  0  THEN  1940 
PRINT  #2,S  :  GOTO  2000 
PRINT  #2 ,  "  1 " ; 

FOR  I  =  O  TO  V  'handles  a  zero  entry  in  result 
PRINT  #2,0; 

NEXT  I 
PRINT  #2,"" 

GOTO  2130 

FOR  I  =  1  TO  NODES 

IF  «(I)  *  0  THEN  2120 

PRINT  #2,  M(I);  'coefficient  of  term 

PTR  *  I 

FOR  J  ■  V  TO  1  STEP  -1  'because  we  are  going 

from  the 

ZZ<J>  =  POT (PTR)  'top  of  the  tree  down 

but  are 

PTR  =  BACK (PTR)  'writing  the  exponents 

in  the 

NEXT  J  'reverse  order 

FOR  J  *  1  TO  V 

PRINT  #2,  ZZ(J); 

NEXT  J 
PRINT  #2,"" 

NEXT  I 
RETURN 

'Addition  subroutine 
PRINT  #2,R0W1;  C0L1;  VI 
FOR  II  *  1  TO  R0W1 
FOR  12  *  1  TO  C0L2 

PRINT  "<";  II;  12;  ">" 

GOSUB  1800  'clear  tree 


'Send  A(I1,I2)  up  through  tree 

FOR  14  *  1  TO  COEFL (II, 12) 

PROD  =  ALPHA  *  A< II  12,14,0) 

FOR  16  -  1  TO  V 

E ( 16)  =  A ( 1 1 , 12 , 14 , 16) 

NEXT  16 
GOSUB  1150 
NEXT  14 

'Send  B(I1,12)  up  through  tree 

FOR  14  =  1  TO  COEFR (II, 12) 

PROD  =  BETA  *  B ( I 1 , 12 , 14 ,0) 

FOR  16  -  1  TO  V 

E  < 16)  =  B ( 1 1 , 12 , 14 , 16) 

NEXT  16 
GOSUB  1150 
NEXT  14 
GOSUB  1870 
NEXT  12 
NEXT  II 
CLOSE  *2 
STOP 

'Matrix  input  subroutine 
INPUT  "Matrix  name";  FILE* 

INPUT  "ROW, COL  dimension";  ROW, COL 
INPUT  "Number  of  polynomial  variates";  V 
OPEN  FILE*  FOR  OUTPUT  AS  #1 

PRINT  "Build  a  square  diagonal  matrix  whose 
diagonal  entries  n; 

PRINT  "are  equal  Y/N";  s  INPUT  OPT1* 

IF  ( (0PT1*  <>  "Y")  AND  (0PT1*  <>  "N"))  THEN  2460 
PRINT  CHR* < 7 ) ;  "***  COEF  EXP1  EXP2  ...  EXPn  ***" 
PRINT  "  a  semicolon  ;  represents  end  of  polynomial 
PRINT  #1,  ROW;  COL;  V 
FOR  I  ■  1  TO  ROW 

FOR  J  ■  1  TO  COL 

S  =  0  :  A*  »  "" 

IF  ( <0PT1*  =  "N">  OR  <<I=1>  AND 
<J=1)))  THEN  2610 

IF  I  =  J  THEN  S  ■  SS  :  PRINT  #1,S  : 

GOTO  2680 
PRINT  #1,1; 

FOR  K  *  0  TO  V  :  PRINT  #1,0;  s  NEXT  K 
PRINT  #1,"H 
GOTO  2750 

PRINT  " <";  STR*(I);  " ,";  STR*(J);  " ) s 
WHILE  (RIGHT* (A*, 1)  <>  ";"> 

S  -  S  +  1 
INPUT  A* 

B* (S)  ■  A*  :  BB* (S)  =  A* 

WEND 

PRINT  #1,S  s  SS  -  S 
FOR  K  =  1  TO  S-l 

IF  0PT1*="Y"  THEN  B*<K)=BB*(K) 
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PRINT  #1 ,BS (K) 

NEXT  K 

IF  0PT1S  =  "Y"  THEN  BS(S)  =  BBS(S> 

PRINT  #1 , LEFTS (BS (S) , LEN (BS (S) ) -1 ) 

PRINT 
NEXT  J 
NEXT  I 
CLOSE  #1 
GOTO  420 

'subroutine  views  a  matrix  file 
INPUT  "Screen  output  Y/N";  0PT1S 

IF  < (0PT1S  <>  "Y")  AND  (0PT1S  <>  "N"))  THEN  2800 

INPUT  "Matrix  to  view";  FILES 

OPEN  FILES  FOR  INPUT  AS  #1 

INPUT  #1,  ROW, COL, V 

IF  OPT IS  =  "Y"  THEN  2890 

LPRINT  "MATRIX  FILES;  ”  IS  “;  STRS(ROW); 

"  X";  STRS (COL) ; 

LPRINT  "  and  has  " ;  STRS(V);  "  polynomial  variates 
GOTO  2920 

PRINT  "MATRIX  " ;  FILES;  "  IS  " ;  STRS (ROW);  "  X"; 
STRS (COL) ; 

PRINT  "  and  has  ";  STRS(V);  "  polynomial  variates" 
PRINT 

FOR  I  ■  1  TO  ROW 

FOR  J  =  1  TO  COL 

IF  OPT IS  *  "Y"  THEN  2970 

LPRINT  FILES;  " (";  STRS(I);  STRS(J);  ”): 

GOTO  2980 

PRINT  FILES;  "(";  STRS(I);  » STRS(J);  "): 
INPUT  #1,S 

FOR  K  *  1  TO  S 

LINE  INPUT  #1,  AS 
IF  OPT  IS  =  "Y"  THEN  3030 
LPRINT  AS  s  GOTO  3040 
PRINT  AS 
NEXT  K 

IF  OPT IS  =  "Y"  THEN  3080 
LPRINT 
GOTO  3090 
PRINT 
NEXT  J 
NEXT  I 
CLOSE  #1 
GOTO  420 

'Subroutine  transposes  a  matrix 

INPUT  "Matrix  to  transpose";  FILES 

INPUT  "Matrix  to  store  transpose  in";  FILEOS 

OPEN  FILES  FOR  INPUT  AS  #1 

INPUT  #1,  ROW, COL, V 

IF  ((ROW  >  ROWS)  OR  (COL  >  COLS))  THEN  1620 
FOR  I  -  1  TO  ROW 


INPUT  #1,  TERMS 
IF  TERMS  >  TERM  THEN  1680 
COEF ( I , J  >  =  TERMS 
FOR  K  =  1  TO  TERMS 
FOR  L  =  O  TO  V 

INPUT  #1,  A ( I , J ,K,L) 

NEXT  L 
NEXT  K 
NEXT  J 
NEXT  I 
CLOSE  #1 

OPEN  FI LEO*  FOR  OUTPUT  AS  #1 
PRINT  #1 ,  COL; ROW; V 
FOR  I  =  1  TO  COL 

FOR  J  =  1  TO  ROW 

PRINT  #1,  COEF (J, I) 

FOR  K  =  1  TO  COEF (J, I) 

FOR  L  =  0  TO  V 

PRINT  #1,  A ( J , I , K , L) ; 

NEXT  L 
PRINT  #1,“" 

NEXT  K 
NEXT  J 
NEXT  I 
CLOSE  #1 
GOTO  420 

'Extracts  U,M,V,N  matrices 

INPUT  "Source  matrix  of  U,M,V,N";  FI LEI* 

INPUT  "Disk  to  output  to  A/B";  D* 

IF  <  (D*  <>  "A;1)  AND  <D*  <>  "B">>  THEN  3490 
D*  -  D*  +  " : " ’ 

PRINT  CHR* < 7 ) ;  "***  WARNING:  files  will  be 
output  to  “ ;  D*; 

PRINT  "U,M,V,N  *** 

INPUT  "OK  Y/N";  OPT* 

IF  OPT *  <>  "Y"  THEN  STOP 
FILE2*  =  D*  +  "U" 

FILE3*  =  D*  +  "M" 

FILE4*  -  D*  +  "V" 

FILES*  =  D*  +  "N" 

OPEN  FI LEI*  FOR  INPUT  AS  #1 
OPEN  FILE2*  FOR  OUTPUT  AS  #2 
OPEN  FILE3*  FOR  OUTPUT  AS  #3 
OPEN  FILE4*  FOR  OUTPUT  AS  #4 
OPEN  FILES*  FOR  OUTPUT  AS  #5 
INPUT  #1, ROW, COL, V 
R0W1  =  ROW/2  :  C0L1  =  COL/2 
IF  ( (R0W1  <>  I NT (R0W1 ) )  OR 
(C0L1  <>  I NT  <C0L1 ) ) )  THEN  1760 

FOR  I  *  2  TO  5  :  PRINT  #1 ,R0W1 ,C0L1 , V  :  NEXT  I 
FOR  I  =  1  TO  ROW 

FOR  J  =  1  TO  COL 
INPUT  #1,S 


FIL  -  5 

IF  (I  <=  R0W1 )  THEN  FIL  =  2 

IF  <  <FIL  =  2)  AND  <J  >  COLD)  THEN  FIL 

IF  ( <FIL  =  5)  AND  (J  <=  COLD)  THEN  FIL 

PRINT  #FIL,S 

FOR  K  =  1  TO  S 

LINE  INPUT  #1,A$ 

PRINT  #FIL,A* 

NEXT  K 
NEXT  J 
NEXT  I 
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That  is,  all  matrices  are  elements  of  the  ring 

cmXn(r  > 

for  m  and  n  of  appropriate  size. 

Because  adding  and  multiplying  matrices  (having  multivariate  polynomial 
entries)  is  tedious  in  practice,  an  Interactive  BASIC  program  is  presented  in  the 
appendix.  This  program,  which  can  be  used  on  a  personal  computer,  permits  the  user 
to  perform  operations  on  matrices  having  multivariate  polynomial  entries.  Via  menu 
selections,  the  user  may  perform 

-weighted  addition  between  two  matrices 
-multiplication  between  two  matrices 

-create  matrices,  with  an  option  of  building  a  diagonal  matrix  whose  diagonal 
entries  are  all  equal 
-view  matrices 
-transpose  a  matrix 

-extract  special  submatrices  (U,M,V,N  of  Chapter  IV)  from  a  given  matrix 
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